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ABSTRACT 

In  the  coupled-dipole  method,  an  arbitrary  particle  is  modeled  as  an 
array  of  N  polarizable  subunits  each  of  which  gives  rise  to  only  electric  dipole 
radiation.  The  total  scattering  is  calculated  by  summing  the  waves  scattered  by 
each  dipolar  subunit  excited  by  the  incident  wave  and  the  waves  of  all  the 
other  dipolar  subunits.  By  accounting  for  the  dipolar  interactions,  the  accuracy 
of  the  scattering  calculations  improves,  but  the  mathematics  become  more 
complicated.  The  matrix  inversion  and  scattering-order  techniques  are  used  to 
solve  for  the  dipolar  interactions. 

The  Clausius-Mosotti  relation  has  been  the  most  widely  used  effective- 
medium  theory  to  relate  the  polarizability  of  the  dipolar  subunits  to  the 
refractive  index  of  the  bulk  particle  that  the  array  represents.  The 
polarizability  of  the  dipolar  subunits  has  been  calculated  using  the  electric 
dipole  coefficient  from  Mie  theory.  This  alternative  expression  for 
polarizability  leads  to  scattering  calculations  that  better  agree  with  Mie  theory. 

Of  all  scattering  angles,  backscattering  is  the  most  sensitive  to  small 
changes  in  particle  size  and  shape.  The  coupled-dipole  method’s  ability  and 
limitations  for  calculating  backscattering  are  demonstrated.  For  particles  with 
size  parameter  less  than  that  associated  with  the  first  backscattering  minimum, 
the  coupled-dipole  method  agrees  favorably  with  Mie  theory.  For  particles 


iv 

with  larger  size  parameters  the  agreement  decreases,  but  accuracy  generally 
improves  by  increasing  the  number  of  dipolar  subunits  in  the  array. 

Backscattering  of  94  GHz  Doppler  radar  by  raindrops  can  be  used  to 
infer  clear  air  velocity;  backscattering  by  ice  crystals  may  provide  similar 
information.  Backscattering  at  94  GHz  by  randomly  oriented  ice  plates  or 
columns  does  not  agree  with  backscattering  by  equal-volume  ice  spheres  for 
size  parameters  greater  than  0.8.  Backscattering  depends  on  zenith  angle  for 
ice  crystals  whose  principle  axes  are  confined  to  the  horizontal  plane.  The 
relationship  between  first  backscattering  minimum  and  size  parameter  varies 
with  particle  shape  and  zenith  angle.  Backscattering  of  vertically  polarized 
light  is  more  sensitive  to  the  presence  of  ice  columns  while  horizontally 
polarized  light  is  more  sensitive  to  ice  plates. 
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Chapter  1 
INTRODUCTION 
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Scattering  and  absorption  of  electromagnetic  waves  are  due  to  the 
heterogeneity  of  matter.  All  matter  is  composed  of  discrete  charged  particles. 
When  an  electromagnetic  wave  excites  a  charged  particle  it  oscillates;  the 
acceleration  of  the  particle  gives  rise  to  electromagnetic  energy  radiated  in  all 
directions.  This  radiation  is  called  scattered  radiation.  If  oscillation  of  the 
particle  is  not  in  phase  with  the  incident  wave,  some  of  the  incoming  energy  is 
absorbed. 

Understanding  this  microscopic  description  of  the  scattering  process  is 
not  a  prerequisite  for  all  scattering  calculations.  For  example,  the  laws  of 
specular  reflection  and  of  refraction  were  known  empirically  long  before  J.  J. 
Thompson  demonstrated  that  atoms  contain  electrons.  In  these  cases  the 
microscopic  actions  of  charged  particles  can  be  described  by  empirically 
derived  macroscopic  equations.  Although  these  macroscopic  solutions  to 
scattering  problems  are  qualitatively  simple  to  understand  and  use,  they  are 
easiest  to  apply  to  simple  geometries  such  as  optically  smooth  planes. 

Scattering  by  macroscopic  particles  can  be  calculated  using  macroscopic 
equations.  If  the  shape  of  a  particle  can  be  described  by  one  of  the  coordinate 
systems  for  which  the  technique  of  separation  of  variables  can  be  exploited,  an 
analytic  solution  for  scattering  may  be  found.  Solutions  exist  for  spheres 
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(Mie,  1908),  infinite  cylinders  (Rayleigh,  1881),  and  spheroids  (Asano  and 
Yamamoto,  1975).  Logan  (1965)  provided  an  interesting  historical  account  of 
the  solution  for  scattering  by  a  sphere. 

The  task  of  calculating  the  scattering  by  an  arbitrary  particle  becomes 
more  difficult;  several  numerical  methods  are  available.  The  T-matrix  method 
is  an  integral  formulation  that  obtains  a  solution  by  iteratively  matching 
boundary  conditions  of  the  internal  and  scattered  fields  of  the  particle.  The 
method  was  first  developed  by  Waterman  (1965)  and  later  made  less  restrictive 
by  Waterman  (1971).  Barber  and  Yeh  (1975)  applied  the  T-matrix  method  to 
scattering  problems  and  renamed  it  the  Extended  Boundary  Condition  Method. 
In  the  point-matching  method  (Oguchi,  1973)  the  fields  inside  and  outside  the 
particle  are  calculated,  and  the  tangential  field  components  are  required  to  be 
continuous  or  matched  across  the  boundary.  Perturbation  methods  treat 
nonspherical  particles  as  spheres  with  a  distorted  boundary  to  estimate  the 
scattered  field  (Yeh,  1964).  The  Purcell-Pennypacker  method  (Purcell  and 
Pennypacker,  1973)  solves  the  problem  of  scattering  by  an  arbitrary  particle  by 
representing  the  particle  as  an  array  of  polarizable  subunits.  The  total 
scattered  field  is  determined  by  summing  the  waves  scattered  by  each  dipolar 
subunit  (also  referred  to  as  dipoles)  that  is  excited  by  the  incident  wave  as  well 
as  by  the  scattered  waves  of  all  the  other  dipoles. 

In  this  dissertation,  the  Purcell-Pennypacker  method  is  used.  Recently, 
it  has  appeared  in  print  under  several  different  names:  the  coupled-dipole 
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approximation  (Singham  and  Salzman,  1986);  the  digitized  Green’s  function 
(Goedecke  and  O’Brien,  1988);  and  the  discrete  dipole  approximation  (Draine, 
1988).  Hereinafter  it  will  be  referred  to  as  the  coupled-dipole  method,  a  name 
that  better  conveys  the  physical  principles  underlying  the  method. 

Gray  (1916),  Kirkwood  (1936),  and  Cassim  and  Taylor  (1965) 
considered  the  interaction  of  dipoles,  i.e.  coupled  dipoles;  however,  Purcell  and 
Pennypacker  were  the  first  to  include  the  interactions  when  calculating 
absorption  and  scattering  by  optically  homogeneous  particles  and  to  compare 
the  results  with  an  analytic  solution.  Having  a  computer  to  perform  the 
mathematical  calculations  was  an  advantage  that  Purcell  and  Pennypacker  had 
over  their  predecessors,  and  the  larger  size  and  faster  speed  of  today’s 
computers  continue  to  attract  scientists  and  engineers  to  the  coupled-dipole 
method.  This  is  evident  when  examining  the  number  of  published  papers 
which  include  calculations  by  or  modifications  to  the  coupled-dipole  method. 

In  the  first  few  years  following  1973,  most  articles  citing  Purcell  and 
Pennypacker  contained  simple  applications  of  the  method  or  only  mentioned 
its  existence.  The  number  of  citations  from  1974  to  1977  totaled  seven;  nearly 
all  were  from  articles  published  in  astrophysical  journals.  Then,  in  an  optics 
journal,  Yung  (1978)  published  an  iterative  technique  for  solving  the  dipole 
interactions.  By  exploiting  the  symmetry  of  a  sphere  and  of  the  incident  plane 
wave,  he  was  able  to  increase  the  number  of  dipolar  subunits  to  over  15,000 
from  the  several  hundred  used  by  Purcell  and  Pennypacker.  Shortly  after  this, 
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Chiappetta  (1980)  published  another  iterative  solution  that  he  referred  to  as  a 
multiple  scattering  technique.  The  advantage  of  this  technique  is  a  substantial 
decrease  of  computer  memory  required  to  carry  out  the  scattering  calculations. 
The  number  of  citations  of  the  original  Purcell-Pennypacker  paper  grew  to  ten 
from  between  1978  and  1980;  between  1981  and  1983  it  decreased  to  three. 

The  popularity  of  them  again  increased  through  the  mid-1980s.  Astronomers 
studying  the  characteristics  of  interstellar  particulate  matter  (Wright,  1988; 
Draine,  1988)  were  joined  by  scientists  interested  in  scattering  by  chiral 
particles  and  dielectric  helices  (e.g.  Singham  et  al.  1986a;  Chiappetta  and 
Torresani,  1988),  agglomerated  soot  particles  (Jones,  1979;  Drolen  and  Tien, 
1987),  dielectric  cubes  (Kattawar  et  al.,  1987),  fractal  dust  grains  (Wright, 

1987),  foreign  objects  on  a  flat  surface  (Taubenblatt,  1990),  and  ice  crystals 
(O’Brien  and  Goedecke,  1988;  Flatau  et  al.  1988,  1990;  Evans  and 
Vivekanandan,  1990;  Vogelmann  et  al.,  1990). 

During  the  last  few  years  several  alternative  solution  techniques  were 
published.  Singham  et  al.  (1986b)  showed  how  the  scattering  properties  of 
randomly  oriented  particles  could  be  calculated  quickly  having  first  determined 
the  inverse  of  the  interaction  matrix.  Singham  and  Bohren  (1988)  reexamined 
the  scattering-order  technique,  which  was  published  earlier  by  Chiappetta 
(1980).  To  reduce  computational  roundoff  error,  Draine  (1988)  added  several 
mathematical  steps  to  the  conjugate  gradient  technique  for  solving  the 
interaction  matrix.  He  also  described  how  the  computation  time  could  be 
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reduced  if  the  particle  shape  with  respect  to  the  incident  plane  wave  was 
symmetric.  Flatau  et  al.  (1990)  reported  on  the  block-Toeplitz  method  for 
inverting  the  interaction  matrix. 


Cloud  and  precipitation  remote  sensing  at  94  GHz  has  become  possible 
only  recently  by  advances  in  microwave  radar  technology  (Lhermitte,  1987). 
The  pulsed  Doppler  radar  enables  the  meteorologist  to  observe  the  motion  and 
growth  of  small  hydrometeors;  this  is  not  possible  with  radars  that  operate  at 
centimeter  wavelengths.  For  certain  rainfall  events,  the  backscattered  signal  is 
used  for  determining  the  clear  air  velocity  (Lhermitte,  1988).  Mie  theory  is 
used  to  calculate  backscattering  by  spherical  raindrops;  however,  it  is  not 
appropriate  for  calculating  backscattering  by  nonspherical  ice  crystals.  The 
coupled-dipole  method  is  used  to  show  that  minima  occur  in  the  backscattered 
signal  from  ice  crystals  as  a  function  of  size  parameter.  As  a  result  of  the 
nonsphericity  of  the  ice  crystals  that  were  modeled,  these  minima  appear  at 
different  size  parameters  for  different  incident  wave  zenith  angles.  The 
minima  also  depend  on  ice  crystal  shape.  Thus,  it  may  be  possible  to  use 
94  GHz  radar  backscattered  signal  to  estimate  the  shape  and  size  distribution 
of  ice  crystals  in  the  1  to  4  mm  size  range.  By  decreasing  the  wavelength  of 
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the  radar,  information  about  smaller  ice  crystals  would  be  available  (Lhermitte, 
1990). 


Organization 

Since  Purcell  and  Pennypacker  first  published  the  coupled-dipole 
method  in  1973,  its  formulation  has  remained  essentially  unchanged.  A 
general  description  of  the  model  is  presented  in  Chapter  2.  Recently,  the 
model  has  been  extended  to  include  particles  of  optically  active  material 
(Singham;  1986)  and  anisotropic  material  (Varadan  et  al.,  1989). 

An  integral  part  of  the  coupled-dipole  method  is  the  effective-medium 
theory  that  provides  the  relationship  between  the  polarizability  of  the  dipolar 
subunits  and  the  refractive  index  of  the  particle  the  array  represents.  The 
most  widely  used  scheme  for  this  is  the  Clausius-Mosotti  relation.  Draine 
(1988)  added  a  radiative  reaction  term  which  guarantees  that  the  polarizability 
will  have  an  imaginary  component  even  if  the  refractive  index  of  the  particle  is 
purely  real.  This  was  done  to  avoid  violation  of  the  optical  theorem.  A 
complex  polarizability  can  also  be  guaranteed  by  using  the  electric  dipole  term 
from  Mie  theory  (Doyle,  1988).  Incorporation  of  this  scheme  in  the  coupled- 
dipole  method  and  consequent  improvements  in  scattering  calculations  are 
described  in  Chapter  3  and  Appendix  A. 
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The  coupled-dipole  method  provides  values  of  the  electric  field 
scattered  by  an  arbitrary  p;irticle.  Chapter  4  describes  how  the  input  data  has 
been  manipulated  to  enhance  the  usefulness  of  the  coupled-dipole  method. 
Derivations  of  equations  for  calculating  observable  scattering  parameters  from 
the  output  data  are  also  presented. 

In  remote  sensing  of  the  atmosphere,  often  only  the  backscattered  signal 
is  readily  available.  Chapter  5  reports  the  advantages,  disadvantages,  and 
limitations  of  using  the  couple-dipole  method  for  determining  backscattering  by 
arbitrary  particles.  In  the  final  section,  application  of  the  coupled-dipole 
method  for  determining  the  backscattered  signal  from  ice  crystals  at  94  GHz 
(3.2  mm)  is  discussed. 
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Chapter  2 

THE  COUPLED-DIPOLE  METHOD 


The  coupled-dipole  method,  which  is  based  on  a  simple  physical 
description  of  how  a  particle  scatters  electromagnetic  waves,  has  remained 
essentially  unchanged  since  it  was  published  by  Purcell  and  Pennypacker 
(1973).  Its  formulation  is  based  on  representing  an  arbitrary  particle  by  an 
array  of  N  dipolar  subunits  arranged  on  a  lattice.  The  dipolar  subunits  are 
sufficiently  small  to  give  rise  to  only  electric  dipole  radiation.  Total  scattering 
is  then  calculated  by  summing  the  electromagnetic  waves  scattered  by  each 
dipolar  subunit  excited  by  the  incident  wave  as  well  as  the  waves  scattered  to 
it  from  all  other  dipolar  subunits. 

The  coupled-dipole  method  can  be  thought  of  as  an  extension  of  the 
Rayleigh-Gans  theory.  Whereas  interactions  among  the  dipolar  subunits  are 
ignored  in  the  Rayleigh-Gans  theory,  they  form  an  integral  part  of  the  coupled- 
dipole  method.  Including  interactions  improves  the  accuracy  of  scattering 
calculations,  but  also  necessitates  solving  a  3N  x  3N  complex  matrix.  Gray 
(1916),  Kirkwood  (1936),  and  Cassim  and  Taylor  (1965)  published  earlier 
investigations  similar  to  the  coupled-dipole  method;  however,  Purcell  and 
Pennypacker  were  apparently  the  first  to  apply  the  method  to  absorption  and 
scattering  by  optically  homogeneous  particles  and  compare  the  results  with  an 
analytic  solution. 
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The  key  to  the  coupled-dipole  method  is  the  interaction  of  the  dipolar 
subunits.  The  following  derivation  of  the  electric  dipole  radiation  field  will 
provide  a  better  understanding  of  these  interactions.  The  derivation  is  similar 
to  that  by  Jackson  (1975,  Chap  9);  however,  these  equations  are  presented  in 
SI  units.  Once  the  interactions  are  known,  the  solution  may  be  obtained  by 
direct  inversion  of  the  3N  x  3N  matrix  or  by  numerical  iteration.  Several 
options  are  discussed  in  this  chapter. 


Formulation  of  the.  Dipole  JnteractiQn.EquatiQ.n 


The  vector  potential  A  is  a  function  of  the  current  density  J  in  the 
system: 


i*0|x-r 


x-x 


3 


dx 


(2.1) 


where  is  the  permeability  of  the  material,  ko  is  the  free  space  wavenumber 
given  by  1^  =  2jt/a0,  a0  is  the  free  space  wavelength,  and  x'  and  "x  are  the 
locations  of  the  charge  and  where  the  field  is  being  evaluated,  respectively.  A 
sinusoidal  time  dependence  exp{-iwt}  of  the  current  is  assumed.  Rapid  spatial 
oscillation  of  the  electromagnetic  field  allows  the  following  approximation  for 
the  argument  of  the  exponential  term 
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|jc - Jc ' |  ~  r-ri'x'  (2-2) 

where  r  =  |  x  j  and  n  is  a  unit  vector  in  the  direction  of  ~x.  For  the  far 
(radiation)  zone  the  distance  from  the  charge  to  the  detector  in  the 
denominator  may  be  replaced  by  r  only.  When  the  source  dimensions  d  are 
small  compared  with  the  wavelength,  exp{-ikon-x'}  may  be  expanded  in 
powers  of  k.  For  electric  dipole  radiation  only  the  first  term  of  this  expansion 
is  kept,  which  simplifies  the  vector  potential  to 

A  -  — — fj(x')dlx'  ,23> 

4 itr  * 

Thus,  the  vector  potential  behaves  like  an  outgoing  spherical  wave  with  a  r'1 
dependence.  In  fact,  equation  2.3  is  valid  everywhere  outside  the  source 
(Jackson,  1975,  Sect  9.2).  The  only  limiting  assumption  so  far  is  d  <?  a0. 

Equation  2.3  can  be  further  manipulated  using  integration  by  parts  and 
by  applying  the  continuity  equation  for  current  density.  The  vector  potential 
becomes 


A  = 


Ho^o  .  eil°r 
4  7 zc  r 


(2.4) 


where  p  is  the  electric  dipole  moment  and  c  is  the  speed  of  light  in  vacuo. 

The  magnetic  induction  is  given  by:  B  =  v  x  A;  while  outside  the  source  the 

— ►  — >  — > 

electric  field  is  E  =  if  v  x  B.  The  electric  field  E  is  now  calculated 

ko 
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E-  J^(Axp)xn^  +  ^^PtIl(±.ih 

4"60  r  4’'e0  ^r3  r2 

Each  dipolar  subunit  has  a  dipole  moment  given  by 

p  -  aE  (2-6) 

where  a  is  the  polarizability  of  the  dipole.  Equations  2.5  and  2.6  are  now 
combined  to  form  a  system  of  3N  equations;  the  resulting  3N  x  3N  matrix  is 
referred  to  as  the  interaction  matrix.  With  N  dipolar  subunits,  the  exciting 
field  E;  at  the  ith  dipole  due  to  the  incident  radiation  and  the  scattered 
radiation  from  the  other  dipolar  subunits  is: 


(2.7) 


where 


eik °r*  [ j2  iko  1 

hr  4tt  z  r  k°+  r  “ 

it7r8ory  y  rl7 


;  b.j-3a.j  2k0> 


(2.8) 


E0  is  the  incident  field,  and  ri}  =  |3fj  -  is  the  distance  from  the  ith  to  the 
jth  dipolar  subunit.  The  polarizability  a  can  be  a  complex  tensor  (Post,  1962), 
but  here  it  is  a  scalar  since  only  isotropic  scatterers  will  be  represented  by 
arrays  of  spherical  dipoles.  The  solution  to  equation  2.7  yields  the  resultant 
electric  field  at  each  dipole  location  with  retardation  taken  fully  into  account. 


The  only  limitation  thus  far  is  that  the  diameter  of  the  dipolar  subunits 


must  be  small  enough  compared  with  the  wavelength.  Since  the  dipolar 
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subunits  are  arranged  on  a  simple  cubic  lattice  and  are  touching  (not 
penetrating)  their  neighbors,  the  dipolar  subunit  diameter  may  be  taken  to  be 
the  same  as  the  lattice  spacing.  Therefore,  the  lattice  spacing  must  also  be 
smaller  than  the  wavelength.  This  is  equivalent  to  the  restriction  required 
when  simplifying  equation  2.3,  i.e.  d  x.  Purcell  and  Pennypacker  (1973) 
expected  good  results  when  kod  <  0.7.  From  the  method-of-moments 
literature  the  grid  spacing  should  not  exceed  one-tenth  of  the  wavelength 
(l^d  <  0.6)  (Massoudi  et  al.,  1984).  Yung  (1978)  placed  a  tighter  restriction  on 
the  spacing,  l^d  <  0.33,  citing  the  necessity  for  smaller  phase  difference  of  the 
incident  wave  between  neighboring  dipolar  subunits.  Recently,  Draine  (1988) 
considered  the  wavelength  within  the  dipolar  subunit  and  included  the  complex 
refractive  index  m  of  the  material.  To  achieve  10%  accuracy  (based  on  a 
zero-frequency  limit)  his  calculations  suggested  kod  j  m  |  <  1.  Table  A.  1 
(page  111)  contains  examples  of  how  the  accuracy  of  scattering  calculations 
using  the  coupled-dipole  method  generally  decreases  as  k^d  |  m  |  increases. 

The  derivation  of  the  coupled-dipole  method  as  described  above  is  a 
heuristic  approach  founded  on  physical  principles.  A  derivation  of  the  method 
based  on  volume  integrals  using  the  free  space  Green’s  functions  has  been 
given  by  Lakhtakia  ( 1990).  The  results  of  both  derivations  are  identical,  and 
the  user  is  left  to  chose  from  several  techniques  to  solve  the  overall  scattering 
problem.  Several  techniques  are  described  in  the  following  section. 
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Solution  Techniques 

Equation  2.7  forms  the  basis  of  the  coupled-dipole  method  by  describing 
the  electric  dipole  field  radiated  by  each  dipolar  subunit.  The  task  is  to  solve 
the  3N  linear  algebraic  equations  to  yield  the  resultant  field  at  all  subunits. 
Several  solution  techniques  are  available,  and  most  are  standard  matrix 
solution  procedures.  Purcell  and  Pennypacker  (1973)  used  an  iterative 
technique  of  successive  substitution;  Yung  (1978)  was  the  first  to  use  the 
conjugate  gradient  technique  for  solving  these  simultaneous  algebraic 
equations;  Singham  et  al.  (1986b)  inverted  the  interaction  matrix;  Chiapetta 
(1980)  and  Singham  and  Bohren  (1988)  used  an  iterative  technique  known  as 
the  scattering-order  technique;  and  Flatau  et  al.  (1990)  exploited  the  block- 
Toeplitz  structure  of  the  interaction  matrix  for  rectangular  particles.  The  two 
solution  techniques  of  interest  here  are  the  matrix  inversion  and  scattering- 
order.  The  advantages  and  disadvantages  of  these  and  some  other  techniques 
are  now  discussed. 


Matrix  Inversion  Technique 

The  system  of  simultaneous  equations  arising  from  equation  2.7  may  be 
expressed  in  the  form  A  =  XB,  where  X  is  the  3N  x  3N  interaction  matrix,  A 
is  the  3N  x  1  vector  containing  the  incident  electric  field  components  E0(r;), 
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and  B  is  the  3N  x  1  vector  containing  the  unknown  exciting  field  components 
E;.  The  matrix  inversion  technique  solves  for  B  by  inverting  the  interaction 
matrix:  X'lA  =  B.  By  the  physical  nature  of  the  problem  an  inverse  exists, 
but  since  the  matrix  inversion  technique  is  a  computer-oriented  routine,  the 
inverse  is  obtained  only  if  the  matrix  is  algorithmically  nonsingular.  As  will  be 
discussed  later,  this  technique  offers  an  advantage  because  iterative  techniques 
rely  on  mathematical  series  that  do  not  always  converge  or  converge  quickly. 
The  inverse  is  obtained  using  UNPACK1  subroutines  CGECO  and  CGEDI. 

Another  advantage  of  the  matrix  inversion  technique  is  that  the 
scattered  field  for  the  particle  in  different  orientations  can  be  calculated 
quickly  (Singham  et  al.,  1986b).  When  the  coordinate  system  describing  the 
particle  remains  fixed,  the  interaction  matrix  and  its  inverse  remain  the  same. 
Thus,  simulation  of  particle  rotation  is  accomplished  by  rotating  the  incident 
field  in  the  opposite  sense. 

The  disadvantage  of  the  matrix  inversion  technique  is  that  it  is  very 
computer  intensive  in  terms  of  storage  requirements  and  computation  time. 

For  a  particle  represented  by  an  array  of  N  dipolar  subunits,  the  3N  x  3N 
complex  interaction  matrix  must  be  stored.  Compiling  the  program  in  single 
precision  (COMPLEX  *  8  on  an  IBM  ES/3090-600)  for  a  250  dipolar  array 
requires  approximately  8.0  megabytes  of  memory.  When  the  program  is 

'LINPACK  is  a  nonproprietary  collection  of  subroutines  amassed  by  the 
National  Bureau  of  Standards  and  available  at  many  computer  facilities. 
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compiled  in  double  precision  rather  than  single  precision,  the  scattering 
calculations  differ  by  less  than  one  percent  for  spheres  with  refractive  indices 
between  1  and  2  and  size  parameters  between  0.5  and  2.0.  Thus,  additional 
memory  burden  is  avoided  by  not  declaring  the  variables  as  COMPLEX  *  16. 
This  is  not  a  general  recommendation  for  not  using  double  precision  because, 
if  the  dipolar  subunits  are  made  very  small  compared  with  the  incident 
wavelength  and  the  incident  field  at  adjacent  dipolar  subunits  are  nearly  equal, 
error  accumulation  due  to  E;  *  Ei+1  may  require  double  precision.  Since 
executable  memory  size  increases  as  N2,  this  matrix  solution  technique  cannot 
be  used  with  small  computers. 

This  solution  technique  may  not  be  acceptable  if  CPU  time  is  a  limiting 
factor.  The  inverse  matrix  is  determined  by  Gaussian  elimination  which, 
although  a  direct  method  of  solution,  is  the  slowest  numerical  method. 
Computer  run  time  to  invert  a  matrix  associated  with  a  160-dipole  array  takes 
approximately  1  minute  on  an  IBM  ES/3090-600,  which  has  a  clock  speed  of 
15  ns  (six  processors  rated  at  about  65  MIPS).  Since  computation  time  is 
proportional  to  N3,  increasing  the  array  to  600  dipoles  will  increase 
computation  time  to  nearly  one  hour.  However,  this  technique  may  be 
attractive  if  the  scattering  response  is  to  be  calculated  for  multiple  orientations 
of  the  particle. 


Scattering-Order  Technique 
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The  second  solution  technique  of  interest  is  the  scattering-order 
technique.  This  approach  was  first  described  by  Chiapetta  (1980)  and 
independently  reexamined  by  Singham  and  Bohren  (1988).  The  advantage  of 
using  the  scattering-order  technique  is  that  the  interaction  matrix  is  not 
required  and  a  much  larger  dipolar  array  can  be  used.  The  scattering-order 
technique  is  based  on  internal  scattering  orders  among  the  dipolar  subunits. 
The  following  is  a  brief  description  of  the  scattering-order  technique. 

For  the  collection  of  N  dipolar  subunits  the  scattering-order  technique 
still  relies  on  equation  2.7,  which  uses  the  incident  field  and  dipole  interaction 
to  calculate  the  electric  field  at  each  dipole;  however,  it  can  be  rewritten  as: 


£i  -  h * E  Cij£j  (2.9) 

i*j 

where  the  matrix  C;j  is  a  function  of  aij,  by,  a ,  and  n^.  Equation  2.9  can  be 
expressed  as  an  infinite  series  in  powers  of  Qj 


**i  j*k 


l+J 


+  EEE<V<V<VV- 

i+j  j+k  k+m 


(2.10) 


Equation  2.10  is  equivalent  to  a  Born  expansion  for  scattering  and  is  described 
as  a  multiple  scattering  solution  in  the  following  fashion.  The  incident  wave 
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hits  the  particle  and  excites  the  dipoles  (Oth  scattering  order  or  the  first  term 
on  the  right  side  of  equation  2.10).  The  dipolar  subunits  then  radiate  to  their 
neighbors  for  the  first  scattering  order  (second  term  of  equation  2.10).  The 
new  field  is  calculated  at  each  dipole  location  and  the  next  scattering  order  is 
calculated--and  so  on  until  the  series  converges.  This  process  can  be  carried 
out  efficiently  by  recognizing  that  each  higher  order  field  may  be  obtained  by 
using  the  results  of  the  previous  scattering  order: 


(2.11) 


i+j  i+j  i*j 

where  Etw  =  £  Cij*Ej(k  l),  k  =  1,  2,  3,  •  •  *.  Hence,  each  scattering  order  can 
«*J 

be  calculated  using  the  previously  determined  fields  and  interaction  matrix. 

The  scattered  field  is  obtained  by  summing  consecutively  higher-order  terms  in 
the  expansion  until  convergence  is  achieved.  In  this  technique  the  field  at  each 
dipole  is  completely  replaced  for  each  scattering  order.  This  reduces  computer 
storage  requirements  since  only  a  few  small  arrays  are  needed  to  calculate 
subsequent  terms  in  the  series;  therefore,  storage  of  the  large  3N  x  3N 
interaction  matrix  is  no  longer  necessary. 

The  disadvantage  of  the  scattering-order  technique  is  that  convergence 
is  not  guaranteed.  If  the  field  amplitude  at  a  dipole  becomes  greater  than  the 
incident  field,  successive  orders  may  continue  to  increase  and  the  series 
ultimately  diverges  (Singham  and  Bohren,  1988).  The  divergent  behavior  of 
the  series  depends  on  the  particle’s  shape,  size,  and  relative  refractive  index. 
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Convergence  is  more  likely  for  elongated  shapes,  smaller  size  parameters,  and 
relative  refractive  indices  that  are  not  much  different  from  unity.  For 
example,  the  series  diverges  for  a  sphere  with  a  relative  refractive  index  of 
1.33  +  iO.O  when  the  size  parameter  exceeds  about  3,  and  for  a  relative 
refractive  index  of  1.9  +  iO.O,  divergence  is  observed  at  a  size  parameter  of 
1.7.  This  places  a  restriction  on  the  particles  that  can  be  studied  using  the 
scattering-order  technique. 

When  the  series  converges  with  little  or  no  oscillation,  the  scattering- 
order  technique  is  relatively  fast.  For  a  160-dipole  array,  one  iteration  takes 
approximately  1  CPU  second  on  an  IBM  ES/3090-600.  Thus,  if  this  series 
converges  in  less  than  60  iterations  the  scattering-order  technique  would  be 
more  attractive  than  matrix  inversion.  By  comparison  the  matrix  inversion 
technique  takes  about  60  seconds  to  solve  a  160-dipole  array.  Table  2.1  lists 
the  number  of  terms  in  the  series  required  to  achieve  convergence  when 
modeling  arrays  of  dipoles  that  represent  a  sphere  with  a  refractive  index  of 
1.9  +  iO.O.  The  convergence  criterion  is  satisfied  when  three  consecutive 
values  of  the  calculated  exciting  field  differ  by  less  than  an  amount  specified 
by  the  user.  A  tolerance  value  of  0.1%  appears  sufficient.  However,  this  does 
not  guarantee  that  the  series  has  converged.  It  is  also  best  to  perform  the 
convergence  test  on  the  backscattered  field  since  convergence  is  slowest  in  the 
backward  direction.  Computation  time  for  the  scattering-order  technique 
varies  as  N2.  Thus,  for  larger  arrays  this  technique  may  be  more  appealing 
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Table  2.1  An  example  of  how  the  number  of  scattering  orders  necessary  for 
series  convergence  increases  with  size  parameter.  The  series  is 
considered  to  have  converged  when  three  successive  values  of 
the  backscattered  field  differ  by  less  than  0.1%. 


Size  Parameter  Number  of  Dipoles  Time  per 

iteration  (sec) 


Number  of 
terms  in  series 


0.7 

461 

9 

7 

1.2 

1064 

18 

19 

1.4 

1357 

23 

24 

1.6 

1791 

30 

50 

20 


than  the  matrix  inversion  technique  except  when  multiple  orientations  of  the 
particle  are  required. 


Other  Techniques 

Two  other  techniques  to  solve  the  interaction  matrix  are  the  conjugate 
gradient  and  the  block-Toeplitz  techniques. 

Yung  (1978)  was  the  first  to  adapt  the  conjugate  gradient  algorithm  for 
the  coupled-dipole  method.  This  is  a  nonlinear  iterative  technique  that  will 
reach  a  solution  in  3N  steps,  but  because  the  algorithm  approaches  the  solution 
quickly,  the  series  can  be  terminated  early  if  a  small  error  is  acceptable 
(Strang,  1986,  Chap  5).  Computer  time  for  an  iteration  is  proportional  to 
N2logN.  Like  the  scattering-order  technique,  the  series  within  the  conjugate 
gradient  technique  will  not  always  converge. 

Draine  (1988)  used  the  conjugate  gradient  technique  to  calculate 
scattering  by  particles  with  fourfold  symmetry.  Exploiting  this  symmetry  will 
reduce  computing  time  by  an  order  of  magnitude  or  more  and  ease  storage 
requirements  by  about  40%.  Unfortunately,  not  all  particles  have  this 
symmetry,  and  even  when  they  do,  this  computational  short-cut  is  limited  to 
only  a  few  particle  orientations. 

Flatau  et  al.  (1990)  introduced  the  block-Toeplitz  technique  into  the 
coupled-dipole  method.  This  technique  exploits  symmetries  of  the  interaction 
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matrix  when  the  dipole  array  represents  a  cube  or  parallelepiped.  Minimal 
storage  requirements  are  the  strong  point  of  the  block-Toeplitz  technique,  and 
while  it  is  much  faster  than  the  matrix  inversion  technique,  it  is  slightly  slower 
than  the  conjugate  gradient  technique. 

All  matrix  solution  techniques  have  advantages  that  make  one  more 
attractive  than  another  depending  on  the  application.  In  this  dissertation  the 
matrix  inversion  and  scattering-order  techniques  were  chosen  for  particular 
purposes.  The  matrix  inversion  technique  is  used  in  Chapter  3  for  testing  a 
new  method  for  determining  the  polarizability  of  the  dipolar  subunits.  It  is 
used  again  in  Chapter  5  for  calculating  scattering  by  particles  in  multiple 
orientations.  The  scattering-order  technique  is  also  used  in  Chapter  5  to 
calculate  scattering  by  particles  that  are  represented  by  arrays  too  large  for  the 
matrix  inversion  technique. 
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Chapter  3 

DOYLE’S  METHOD  FOR  CALCULATION  OF  POLARIZABILITY 


Purcell  and  Pennypacker  (1973)  formulated  the  coupled-dipole  method 
to  study  scattering  by  arbitrary  particles  at  non-zero  frequencies.  The  arbitrary 
particle  is  represented  by  an  array  of  polarizable  subunits  arranged  on  a  simple 
cubic  lattice.  To  solve  the  scattering  problem  one  must  choose  a  technique  to 
solve  the  system  of  linear  equations  and  also  determine  how  to  represent  the 
polarizable  subunits.  Several  solution  techniques  were  discussed  in  Chapter  2; 
in  this  chapter  a  scheme  to  determine  the  electric-dipole  polarizability  of  the 
dipolar  subunits  is  described. 

Splitting  the  homogeneous  scatterer  into  smaller  objects  alters  its 
electromagnetic  description.  To  relate  the  discrete  and  continuous  aspects  of 
the  particle,  Purcell  and  Pennypacker  used  the  Clausius-Mosotti  (CM)  relation. 
Like  every  effective-medium  theory  the  CM  relation  is  approximate,  but  it  has 
been  used  by  nearly  all  adherents  of  the  coupled-dipole  method.  Here,  based 
on  work  by  Doyle  (1989),  an  alternative  method  for  calculating  polarizability  is 
presented. 

No  material  is  homogeneous  in  the  absolute  sense.  A  medium  which  is 
perceived  as  homogeneous  is  actually  comprised  of  heterogeneities,  i.e.  charged 
particles.  As  long  as  the  heterogeneities  (referred  to  hereinafter  as  grains)  are 
much  smaller  than  the  wavelength,  the  behavior  of  an  electromagnetic  wave 
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passing  through  such  a  medium  can  be  described  statistically  by  estimating  the 
bulk  optical  property  of  the  particle  using  an  effective-medium  theory  (Bohren 
and  Wickramasinghe,  1977).  An  effective-medium  theory  requires  knowledge 
of  the  volume  fraction  of  the  embedded  grains  and  the  optical  properties  of 
the  grains  and  the  surrounding  matrix. 

When  the  grains  are  composed  of  many  molecules,  the  material  can  be 
characterized  by  an  effective  dielectric  function  using  the  Maxwell  Garnett 
(MG)  theory  (Maxwell  Garnett,  1904;  Bohren  and  Huffman,  1983,  Sect  8.5): 


(l-/)e„+/Pe 

1-/+/P 


(3.1) 


where  eav,  e,  and  era  are  the  dielectric  functions  of  the  bulk  medium,  grains, 
and  surrounding  matrix,  6  is  a  function  of  grain  shape,  and  f  is  the  volume 
filling  factor  of  the  grains.  If  the  grains  are  spherical,  B  =  3em/(e+em),  and 
equation  3.1  can  be  rewritten  as 
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and  the  dielectric  function  of  the  grains  e  must  be  determined.  Rearrangement 
of  equation  3.2  yields 


£  .  2(1 -/)«*„ -Q +/)«„«,  (3 3 

(l-/)em-(l+2/)c„  ' 

This  form  of  the  MG  theory  will  be  mentioned  later  in  the  chapter. 

When  the  grains  are  too  small  to  be  assigned  a  dielectric  function,  the 
Clausius-Mosotti  (CM)  relation  (also  known  as  the  Lorentz-Lorenz  equation) 
may  be  used  to  determine  the  polarizability  a  of  the  grains 


g  .  _  3(m*-l)  (34) 

N(ew  +  2e.)  N(/m2  +  2) 

where  N  is  the  number  density  of  the  grains  per  unit  volume  (N  «  f'1),  and  m  is 
the  complex  relative  refractive  index.  Even  though  the  MG  theory  and  CM 
relation  give  different  physical  properties,  they  are  very  similar.  Barker  ( 1973) 
published  a  derivation  of  the  MG  equation  that  is  formally  equivalent  to  the 
derivation  of  the  CM  relation. 

In  the  coupled-dipole  method  the  grains  are  the  dipolar  subunits.  The 
means  to  estimate  the  dielectric  function  or  the  polarizability  of  the  subunits 
are  given  in  equations  3.3  and  3.4.  The  relationship  between  e  and  a  of  the 
dipolar  subunits  is  obtained  by  combining  these  equations  to  produce  another 
form  of  the  Clausius-Mosotti  relation: 
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^0  -  g  ) 

a  -  4ua3- - (3.5) 

(«+2ew) 

where  a  is  radius  of  the  dipolar  subunit.  This  brief  discussion  of  effective- 
medium  theory  provides  a  basis  for  introducing  Doyle’s  method  for  calculating 
the  electric  dipole  polarizability  of  the  dipolar  subunits. 

For  a  sphere  in  the  longwave  limit  (x<  1,  |  m  |  x  «  1,  where  x  is  the  size 
parameter:  x  =  2jra/A0,  with  a  being  the  radius  and  A0  the  wavelength  of 
electromagnetic  radiation  in  free  space)  the  scattering  amplitude  St  is  given  by 
(Bohren  and  Huffman,  1983,  Sect  5.2) 


(3.6) 


Since  the  sphere  is  in  the  Rayleigh  regime  for  scattering,  the  scattering 
amplitude  can  be  expressed  as  S1  =  (3/2)a1,  where  at  is  the  electric  dipole 
coefficient  from  Mie  theory.  Thus,  the  polarizability  can  be  written  in  terms  of 
the  electric  dipole  coefficient 


a 


(3.7) 


Since  interactions  between  the  dipolar  subunits  in  the  coupled-dipole  method 
are  due  to  electric  dipole  radiation,  this  expression  may  be  appropriate  for 
calculating  the  polarizability  of  the  subunits.  Equation  3.7  does  not  differ  from 
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equation  3.5  in  the  zero  frequency  limit.  If  at  is  expanded  in  terms  of  the  size 
parameter  x  and  the  series  truncated  after  the  first  term,  then  equation  3.7  is 
identical  to  equation  3.5. 

Doyle  (1989)  used  equation  3.7  for  determining  the  reflectance  of  a 
composite  of  silver  spheres  in  a  transparent  medium  and  obtained  calculated 
values  that  agreed  better  with  experimental  values  than  if  only  the  first  term  in 
the  series  had  been  used.  The  task  now  is  to  use  this  expression  in  the 
coupled-dipole  method  and  see  if  it  improves  scattering  calculations. 

From  this  point  on  the  calculation  of  the  electric  dipole  polarizability 
using  will  be  known  as  Doyle’s  method;  equation  3.7  will  be  referred  to  as 
Doyle’s  expression.  The  value  for  a1?  the  electric  dipole  coefficient,  is  given 


-  ^(x^C/nx) 
al  -  - ; - - - 

^(x)ij r^mx) 

where  V>  and  £  are  the  Ricatti-Bessel  functions  given  by: 


(3.8) 


i|f  L(x)  -  — —  -  cos(x)  ,  ijr^x) - +  ------  +  sin(x) 


(3.9) 


$j(x)  -  -i- — ,  ^(x) 
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x 
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Thus,  Doyle’s  expression  gives  the  exact  value  for  the  electric  dipole 
polarizability. 

Most  applications  of  the  coupled-dipole  method  use  the  CM  relation  as 
it  appears  in  equation  3.4.  In  this  form  only  the  dielectric  functions  of  the  bulk 
particle  and  the  surrounding  medium  are  required.  To  use  Doyle’s  method, 
the  complex  refractive  index  m  of  the  dipolar  subunits  must  be  determined  for 
equation  3.8  using  equation  3.3  where  m=(e/em)M. 

Scattering  calculations  using  Doyle’s  expression  to  calculate 
polarizability  in  the  coupled-dipole  method  are  discussed  in  Appendix  A 
beginning  on  page  104.  The  results  are  compared  with  two  other  methods  for 
calculating  the  polarizabilities:  the  CM  relation  and  an  expression  published 
by  Draine  (1988).  Draine  introduced  a  radiative  reaction  term  to  the  CM 
relation  which  guarantees  an  a  non-zero  imaginary  part  of  the  polarizability  to 
account  for  attenuation  when  the  refractive  index  of  the  material  is  real.  The 
solutions  using  these  three  methods  were  compared  with  Mie  theory  (Bohren 
and  Huffman,  1983,  Appendix  A).  In  all  cases  scattering  calculations  using 
Doyle’s  expression  agreed  with  Mie  theory  better  than  when  using  the  CM 
relation  or  Draine’s  radiative  reaction  term. 

Another  test  for  Doyle’s  expression  occurs  at  a  Frohlich  frequency. 

This  is  an  absorption  mode  where  the  relative  complex  dielectric  function 
approaches  -2:  e/  em  -*■  -2  +  i0.  At  the  Frohlich  frequency  the  CM  relation 
yields  an  unbounded  polarizability  while  Doyle’s  expression  remains  bounded. 


A  Frohlich  frequency  occurs  for  quartz  in  the  infrared  part  of  the 
electromagnetic  spectrum.  Scattering  calculations  using  the  coupled-dipole 
method  with  Doyle’s  method  were  shown  to  agree  favorably  with  measured 
data.  A  detailed  discussion  and  results  of  using  Doyle’s  expression  at  a 
Frohlich  frequency  is  found  in  Appendix  A  beginning  on  page  112. 


Chapter  4 

SCATTERING  CALCULATIONS  USING 
THE  COUPLED-DIPOLE  METHOD 

The  description  of  the  coupled-dipole  method  in  Chapters  2  and  3  has 
been  focused  on  its  formulation,  solution  techniques  for  determining  dipolar 
interactions,  and  a  new  scheme  for  determining  the  polarizability  of  the  dipolar 
subunits.  This  chapter  contains  information  on  how  the  input  data  can  be 
manipulated  to  enhance  the  usefulness  of  the  program  and  includes  derivations 
of  some  of  the  observable  scattering  parameters  that  are  calculated  by  the 
program.  A  version  of  the  coupled-dipole  method  computer  program  was 
obtained  from  Singham,  and  the  modifications  described  in  this  chapter  were 
made  to  serve  the  needs  of  the  research  presented  here  and  in  the  future. 


An  advantage  of  the  coupled-dipole  method  over  other  methods  used  to 
calculate  scattering  by  nonspherical  particles  is  that  a  particle  of  any  shape  can 
be  modeled.  In  the  coupled-dipole  method  a  particle  is  represented  by  an 
array  of  subunits  located  on  a  cubic  lattice.  However,  caution  is  required 
when  constructing  the  array  or  choosing  the  coarseness  of  the  lattice  spacing. 
For  example,  consider  an  array  of  eight  dipolar  subunits  (hereinafter  referred 
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to  as  dipoles)  situated  on  the  vertices  of  a  cube  centered  about  the  origin. 

This  simple  array  represents  a  cube;  however,  as  will  be  described  soon,  it 
could  also  represent  a  sphere  constructed  of  all  the  dipoles  within  /3/2  units  of 
the  origin.  This  ambiguity  of  shape  must  be  minimized  when  building  arrays 
of  dipoles.  Another  example  of  poor  particle  representation  is  the  depiction  of 
a  cylinder  as  a  single  string  of  dipoles.  This  configuration  lacks  adequate  cross 
sectional  area  and  no  sagittal  dipole  interaction  necessary  to  represent  the 
cross  section.  Subroutines  for  several  particle  shapes  are  available  in  the 
coupled-dipole  method;  they  are  briefly  described  here,  and  their  listings  are 
found  in  Appendix  B. 

The  unit  of  length  used  to  describe  the  particle  dimensions  in  this 
dissertation  is  the  lattice  space  (Is).  One  lattice  space  is  the  distance  from  the 
center  of  one  dipole  to  the  center  of  an  adjacent  dipole.  For  determining  the 
equivalent  volume  of  a  dipolar  array,  one  dipole  is  considered  to  have  unit 
volume. 

Subroutine  SPHERE  builds  a  dipolar  array  to  represent  a  spherical 
particle  in  the  following  manner.  To  model  a  sphere  of  radius  a,  the 
subroutine  builds  a  cube  with  sides  2a,  centered  at  the  origin.  All  lattice  points 
within  the  cube  that  fall  within  a  distance  of  the  origin  are  included  in  the 
array.  The  position  of  the  origin  with  respect  to  the  lattice  is  optional; 
however,  two  logical  choices  exist.  The  origin  can  be  centered  on  one  dipole 
which  will  be  referred  to  as  the  central  dipole  (0.0,  0.0,  0.0),  or  the  origin  can 
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lie  in  the  midst  of  a  core  of  eight  dipoles  that  are  offset  by  one-half  dipole 
diameter  into  respective  octants  of  the  coordinate  system  (±0.5,  ±0.5,  ±0.5). 
Many  published  papers  use  only  the  latter  option,  which  reduces  the  number  of 
arrays  that  can  be  used  to  represent  a  sphere.  A  sphere  of  each  type  is  shown 
in  Figure  4.1.  The  top  array  represents  a  sphere  with  the  origin  centered  on  a 
dipole  and  contains  251  dipoles.  It  is  constructed  by  choosing  a  radius  of  3.8  Is 
and  has  an  effective  radius  of  3.91  Is.  The  bottom  sphere  has  the  origin 
centered  in  a  core  of  eight  dipoles  and  contains  304  dipoles.  It  is  constructed 
by  choosing  a  radius  of  4.1  Is  and  has  an  effective  radius  of  4.17  Is.  More  will 
be  said  about  the  effective  radius  of  particles  in  the  section  describing  particle 
size,  which  begins  on  page  38.  Scattering  by  a  sphere  calculated  by  the 
coupled-dipole  method  has  been  shown  to  compare  favorably  with  Mie  theory 
(e.g.,  Yung,  1978;  Singham  and  Bohren,  1988). 

Subroutine  SPHEROID  is  similar  to  subroutine  SPHERE  except  that 
criteria  for  including  dipoles  in  the  array  are  based  on  the  principal  radii  of  a 
spheroid.  The  subroutine  is  set  up  to  build  a  prolate  (oblate)  spheroid  with  its 
major  (minor)  axis  in  the  z-direction.  The  spheroid  is  automatically  centered 
around  the  origin;  thus  the  central  dipole  may  or  may  not  be  located  at  the 
origin  depending  on  whether  the  major  and  minor  radii  are  represented  by  an 
odd  or  even  number  of  dipoles.  For  this  reason  subroutine  SPHERE  is  a 
separate  entity  and  not  a  special  case  within  subroutine  SPHEROID. 

Examples  of  spheroids  are  shown  in  Figure  4.2.  The  top  array  represents  a 
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Figure  4.1 


Arrays  of  dipoles  arranged  to  represent  spheres.  The  top  array 
contains  251  dipoles  and  has  an  effective  radius  of  3.91  Is;  the 
bottom  array  contains  304  dipoles  and  has  an  effective  radius  of 
4.17  Is. 
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Figure  4.2  Arrays  of  dipoles  arranged  to  represent  spheroids.  The  top  array 
contains  290  dipoles  and  represents  a  2x  1  oblate  spheroid;  the 
bottom  array  contains  128  dipoles  and  represents  a  2x  1  prolate 
spheroid. 
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2x1  oblate  spheroid,  contains  290  dipoles,  and  has  an  equivalent  spherical 
radius  of  4.11  Is.  The  bottom  array  represents  a  2x  1  prolate  spheroid, 
contains  128  dipoles,  and  has  an  equivalent  spherical  radius  of  3.21  Is.  Unlike 
a  sphere,  scattering  by  spheroids  and  other  particle  shapes  that  will  be 
discussed  is  orientationally  dependent.  Scattering  by  spheroids  using  the 
coupled-dipole  method  has  been  shown  to  compare  favorably  with  the  T- 
matrix  method  at  low  frequencies  and  low  refractive  indices  (Goedecke  and 
O’Brien,  1988). 

The  subroutine  RECTSLD  builds  a  dipolar  array  that  represents  a 
rectangular  solid.  The  user  specifies  the  x-,  y~,  and  z-dimensions,  and  the 
>  subroutine  centers  the  particle  around  the  origin.  This  subroutine  was  written 
to  compare  results  published  by  Purcell  and  Pennypacker  (1973).  The  block- 
Toeplitz  technique  for  solving  the  interaction  matrix  equation  specifically 
handles  particles  of  this  shape  (Flatau  et  al.,  1990).  Unfortunately,  few 
applications  in  meteorology  require  knowledge  of  scattering  by  rectangular 
particles.  An  example  of  this  type  of  particle  is  shown  in  the  top  part  of 
Figure  4.3.  The  array  is  constructed  with  x-,  y-,  and  z-dimensions  of  4,  4,  and 
10,  and  contains  160  dipoles. 

Subroutine  CYLNDER  constructs  a  cylinder  with  its  longitudinal  axis  in 
the  x  direction.  Input  parameters  for  assembling  the  array  include  the  length, 
the  radius  of  the  circular  cross  section,  and  whether  the  center  dipole  of  the 
cross  section  is  located  on  the  x-axis  (x,  0.0,  0.0)  or  if  four  dipoles  are 
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Figure  4.3  Arrays  of  dipoles  arranged  to  represent  a  rectangular  solid  and  a 
cylinder.  The  top  array  contains  160  dipoles  and  represents  a 
4  x  4  x  10  rectangular  solid;  the  bottom  array  contains  240  dipoles 
and  represents  a  cylinder  with  a  cross-sectional  radius  of  2.9  Is 
and  a  length  of  10  Is. 
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centered  around  the  x-axis  [(x,  ±0.5,  ±0.5)  similar  to  the  offset  option  in 
SPHERE].  An  example  of  a  cylindrical  particle  is  shown  in  the  bottom  half  of 
Figure  4.3.  The  array  is  constructed  by  choosing  a  cross-sectional  radius  of 
2.9  Is,  has  a  length  of  10  Is,  contains  240  dipoles,  and  has  an  effective  circular 
cross-section  radius  of  3.04  Is.  Scattering  by  a  finite  cylinder  was  compared 
with  the  results  of  an  analytic  solution  for  infinite  cylinders,  BHCYL  (Bohren 
and  Huffman,  1983,  Appendix  C).  Approximations  based  on  Fraunhofer 
diffraction  theory  indicate  that  cylinders  with  an  aspect  ratio  greater  than  10:1 
may  be  regarded  as  effectively  infinite  (Bohren  and  Huffman,  1983, 

Sect  8.4.7).  For  such  a  cylinder,  scattering  calculations  using  the  coupled- 
dipole  method  and  BHCYL  are  in  good  agreement.  For  example,  values  of 
extinction  efficiency  per  unit  length  for  cylinders  with  a  refractive  index  of 
1.33  +  i0.4  are  within  2%. 

Subroutine  FCC  builds  a  sphere  similar  to  SPHERE  except  that  the 
dipoles  are  placed  on  a  face-centered  cubic  lattice  rather  than  simple  cubic. 
The  creation  of  FCC  was  to  test  a  different  dipole  spacing  scheme.  Since  a 
face-centered  cubic  lattice  has  12  nearest  neighbors  rather  than  only  6  as  with 
the  simple  cubic,  and  it  has  a  packing  fraction  higher  by  a  factor  of  J2  (Kittel, 
1976),  it  was  thought  that  a  better  representation  of  a  sphere  could  be 
achieved  with  little  effort.  The  disadvantage  of  this  structure  would  be  that 
with  an  increased  number  of  dipoles,  more  computer  storage  and  computing 
time  is  required  to  solve  the  interaction  matrix  for  an  equal-sized  sphere. 
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Another  concern  with  this  structure  is  that  mathematical  approximations  that 

hold  for  the  simple  cubic  lattice  may  not  be  valid  for  the  face-centered  cubic 

lattice.  Lorentz  showed  that  for  atoms  on  a  simple  cubic  lattice  the  electric 

field  resulting  from  the  contributions  of  the  nearest  dipoles  (E^)  vanishes  at 

any  lattice  site  (Jackson,  1975,  Section  4.5).  Although  face-centered  cubic  is 

not  equivalent  to  simple  cubic,  the  structure  may  be  symmetric  enough  to 

assume  «  0.  If  this  approximation  were  not  valid  the  CM  relation  (and 

Doyle’s  expression)  would  need  to  be  reformulated.  Scattering  calculations 

were  made  using  SPHERE  and  FCC  and  were  compared  with  results  from 

Mie  theory.  Since  SPHERE  compared  more  favorably  with  Mie  than  FCC  it 

— ► 

appears  that  the  approximation  «  0  may  not  be  valid. 

Subroutine  HEXAGON  constructs  a  hexagon  crystal  with  its 
longitudinal  axis  lying  in  the  x  direction  and  two  vertices  aligned  with  the  z 
axis.  Depending  on  the  input  parameters,  the  array  can  represent  either  a 
hexagonal  plate  or  column.  Because  the  array  is  positioned  on  a  simple  cubic 
lattice,  a  regular  hexagon  is  difficult  to  construct.  Instead  of  120°  angles,  the 
hexagon  is  constructed  with  90°  and  135°  angles,  and  the  cross-section 
dimensions  are  chosen  to  best  represent  an  equilateral  polygon.  Figure  4.4 
contains  two  examples  of  hexagonal  particles.  The  top  figure  is  a  column  and 
has  an  effective  cross-sectional  diameter  of  5.98  Is,  a  length  of  12  Is  or  an 
aspect  ratio  of  about  2x1,  and  contains  336  dipoles.  The  bottom  figure  is  a 


Figure  4.4  Arrays  of  dipoles  arranged  to  represent  hexagonal  crystals.  The 
top  array  contains  336  dipoles  and  represents  a  2x  l  hexagonal 
column,  the  bottom  array  contains  212  dipoles  and  represents  a 
2x  l  hexagonal  plate. 
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plate  and  has  an  effective  cross-sectional  diameter  of  8.21  Is,  a  depth  of  4  Is  or 
an  aspect  ratio  of  about  2x1,  and  contains  212  dipoles. 

Particle  .Size 

The  subroutines  described  in  the  previous  section  can  be  used  to 
construct  arrays  of  dipoles  that  represent  particles  of  various  shapes. 

However,  as  a  result  of  the  discrete  nature  of  the  dipoles,  only  a  discrete  set  of 
particle  sizes  can  be  represented  unless  additional  modifications  are  made. 

The  modifications  to  vary  the  particle  size  are  based  on  either  the  effective 
radius  or  volume  of  particle  that  the  array  represents.  If  a  particle  of  given 
size  is  to  be  modeled,  the  subroutines  adjust  the  radii  of  the  dipoles 
accordingly.  The  following  exemplifies  these  modifications. 

If  SPHERE  were  used  to  build  a  particle  of  radius  3.0  Is  with  the  origin 
at  a  central  dipole,  the  resulting  dipolar  array  would  contain  all  the  dipoles 
within  3.0  Is  of  the  origin  and  consist  of  123  dipoles.  The  volume  of  the 
scatterer  is  then  N,  the  number  of  dipoles.  The  effective  radius  of  the  particle 
is  given  by 
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(4.1) 


In  this  case  ae  =  3.09  Is.  If  the  input  parameter  for  radius  had  been  3.1  Is,  the 
subroutine  would  have  built  the  same  123-dipolar  array.  The  next  largest  array 
would  not  occur  until  the  input  radius  reached  3.2  Is  which  would  result  in  an 
array  of  147  dipoles;  the  effective  radius  for  this  sphere  is  3.27  Is.  Thus,  a 
particle  of  radius  between  3.09  and  3.27  Is  cannot  be  built  without  further 
modifications  to  the  subroutine.  One  modification  has  already  been  discussed: 
the  option  for  repositioning  the  origin  with  respect  to  the  lattice  points.  This 
option  approximately  doubles  the  number  of  arrays  that  can  represent  a 
sphere,  but  overall  selection  remains  limited.  Table  4.1  contains  the  complete 
list  of  spheres  (of  radius  3.0  <  a  <  4.5)  that  can  be  represented  using  the 
SPHERE  subroutine.  The  list  contains  information  on  arrays  that  have  the 
origin  at  a  central  dipole  or  within  a  group  of  eight  dipoles.  It  is  observed  in 
Table  4.1  that  a  spherical  particle  with  an  effective  radius  of  3.25  Is  cannot  be 
built.  The  following  modification  eliminates  this  restriction. 

When  a  particle  of  specific  size  is  required,  the  user  is  required  to  input 
the  actual  dimensions  and  the  radius  of  the  dipole  is  adjusted  so  the  effective 
radius  matches  the  desired  value.  For  instance,  if  a  sphere  with  an  effective 
radius  of  3.25  Is  is  to  be  modelled,  the  sphere  containing  147  dipoles  may  be 
used,  but  the  dipole  radius  will  be  reduced  to  0.48  Is  from  0.50  Is.  The  size  of 
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Table  4.1  Physical  characteristics  of  dipolar  arrays  that  are 

constructed  by  subroutine  SPHERE.  Radius  is  the  input 
parameter,  eff  rad  is  the  effective  radius  based  on 
equivalent  volume  sphere,  and  N  is  the  number  of  dipoles 
in  the  array. 


Origin  not  at  a  central  dipole 


Radius  Eff  Rad  N 


3.0 

3.19 

136 

3.3 

3.37 

160 

3.6 

3.68 

208 

3.9 

4.06 

280 

4.1 

4.17 

304 

4.4 

4.41 

360 

Origin  at  a  central  dipole 


Radius  Eff  Rad  N 


3.0 

3.09 

123 

3.2 

3.27 

147 

3.5 

3.50 

179 

3.8 

3.91 

251 

4.2 

4.18 

305 

4.3 

4.33 

341 
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a  spheroid  is  adjusted  based  on  the  radius  of  an  equivalent  sphere 
(equation  4.1).  For  a  cylinder  and  hexagon  the  modification  depends  on  the 
effective  cross-sectional  radius: 


a 


e 


(4.2) 


In  this  manner  particles  of  any  size  can  be  represented. 

This  ability  to  change  the  size  of  the  dipoles  is  also  important  because 
backscattered  radiation  is  highly  dependent  on  particle  shape.  When  more 
dipoles  are  used  to  represent  a  particle,  the  array  becomes  a  truer 
representation  of  that  particle.  As  will  be  discussed  in  Chapter  5, 
backscattering  calculations  can  be  improved  by  modelling  particles  with  larger 
arrays  of  smaller  dipoles. 


Scattering  Matrices 

Two  related  matrices  are  encountered  when  calculating  scattering  of 
polarized  light  by  a  particle:  the  amplitude  scattering  matrix  and  the  Mueller 
matrix.  The  amplitude  scattering  matrix  provides  the  relation  between  the 
incident  and  scattered  electric  fields.  The  Mueller  matrix  is  an  expansion  of 
the  amplitude  scattering  matrix  using  a  set  of  observable  properties  of 
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polarized  light  which  are  referred  to  as  the  Stokes  vectors.  Both  matrices  are 
calculated  in  the  coupled-dipole  method. 

From  this  point  on  light  will  refer  to  plane  wave  radiation  propagating 
in  the  +z  direction.  Polarization  in  the  context  of  polarized  light  is  not  to  be 
confused  with  the  polarizability  of  the  dipolar  subunits.  The  state  of 
polarization  describes  the  behavior  of  the  electric  field  in  the  jt-y  plane  as  the 
wave  propagates  in  the  z  direction.  For  example,  parallel  linearly  polarized 
light  refers  to  a  linearly  polarized  plane  wave  whose  electric  field  remains 
parallel  to  some  arbitrary  axis.  Conventions  for  defining  polarized  light  are 
not  standard;  the  convention  adopted  here  is  that  used  by  Bohren  and 
Huffman  (1983,  Chap  2.11). 

Singham  et  al.  (1986a)  analyzed  the  scattered  light  in  terms  of  left-  and 
right-circularly  polarized  light.  The  equations  necessary  to  evaluate  parallel 
and  perpendicular  linearly  polarized  light  have  been  added  to  the  coupled- 
dipole  method;  the  derivation  of  these  equations  used  to  calculate  the 
amplitude  scattering  matrix  is  now  described. 

The  amplitude  scattering  matrix  relates  the  parallel  and  perpendicular 
components  of  the  incident  (i)  and  scattered  (s)  electric  fields 
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where  Sj(j  =  1,4)  is  the  amplitude  scattering  matrix  (Bohren  and  Huffman,  1983, 
Sect  3.3).  By  expressing  the  parallel  and  perpendicular  components  of  the 

.  .  — >  — >  — >  — >  — *  /y, 

incident  field  in  terms  of  x  and  y:  E(i  =  Ex,  and  Exl  =  -Ey,  where  Ex  =  EqX» 

Ey  =  E0y’  anc*  E0  =  1,  the  x  and  y  component  values  for  the  incident  field 
become:  Ej-  =  -Ey  =  1  and  EJj  =  Exi  =  0.  With  this  information,  the 
matrix  elements  in  equation  4.3  can  be  solved  for  in  terms  of  x  and  y 
components  of  parallel  and  perpendicular  linearly  polarized  scattered  light: 


where  C  =  exp{ikr}/-ikr. 

The  electric  field  at  a  detector  Ed  placed  at  distance  rd  in  the  unit 
direction  is  obtained  by  summing  the  far-field  amplitudes  from  the  N 
dipoles 


u2  *Td  _  5 

(1  -  ndnd)  £  e  Ej . 

Td  7-1 


N 


(4.5) 


where  1  is  the  identity  matrix.  Equation  4.5  is  similar  to  equation  4.3.  With 
minor  algebraic  manipulation  the  elements  of  the  amplitude  scattering  matrix 
can  be  calculated  using  the  coupled-dipole  method. 
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The  values  for  the  amplitude  matrix  are  passed  into  a  subroutine  where 
the  4x4  Mueller  matrix  is  defined  by 


K  sa  sa  sl4||  i, 

Qs  ^21  ^22  $23  $24  Qt  ^  ^ 

us  "  s3l  S32  S33  su  ut  * 

Vs)  1^41  *42  *43  *JU, 

(Bohren  and  Huffman,  1983,  p  65).  The  Mueller  matrix  describes  the 
relationship  between  the  incident  and  scattered  Stokes  vector.  The  Stokes 
vector  (I,  Q,  U,  and  V)  is  a  description  of  the  state  of  polarization  of  both 
incident  (i)  and  scattered  (s)  light.  Monochromatic  light  of  any  state  of 
polarization  can  be  represented  by  the  Stokes  vector.  The  I  component 
represents  the  irradiance  of  the  wave;  V  defines  its  handedness.  Q  is  related 
to  horizontally  and  vertically  polarized  light,  and  0  is  related  to  ±45°  linearly 
polarized  light.  Horizontal  and  vertical  are  defined  with  respect  to  the 
direction  of  wave  prop°gation.  (See  Bohren  and  Huffman,  1983, 

Section  2.11.1,  for  a  more  detailed  derivation  of  the  Stokes  vectors.) 

In  Chapter  5,  the  scattering  and  detection  of  horizontally  and  vertically 
polarized  light  will  be  discussed.  When  a  plane  wave  has  been  linearly 
polarized  by  an  ideal  linear  polarizer,  the  electric  field  vector  remains  parallel 
to  a  particular  axis  called  the  transmission  axis.  The  Mueller  matrix  for  an 
ideal  linear  polarizer  is 
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2 


1  cos  2^  sin2£  0^ 

cos25  cos22^  cos2£  sin25  0 
sin25  sin25  cos2£  sin22£  0 
0  0  0  0 


(4.7) 


where  £  is  the  angle  between  the  transmission  axis  and  the  horizontal  plane 


(horizontal  with  respect  to  the  direction  of  wave  propagation).  When  a 
particle  is  illuminated  with  horizontally  polarized  light,  the  amount  of 
backscattered  light  that  is  horizontally  polarized  can  be  determined  by 
multiplying  a  series  of  Mueller  matrices;  these  are  shown  in  equation  4.8. 
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^10  0^ 
110  0 
0  0  0  0 
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(1  1  0  0^ 
110  0 
0  0  0  0 
0  0  0  0 


fl' 

0 

0 
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(4.8) 


The  4x1  matrix  on  the  right  is  the  incident  Stokes  vector  representing 
normalized  irradiance.  Next  to  the  Stokes  vector  is  the  Mueller  matrix  for  a 
linear  polarizer  with  a  horizontal  transmission  axis  (£  =  0°).  The  product  of 
these  two  matrices  yields  the  4  x  1  Stokes  vector  that  represents  horizontally 
polarized  light.  (Sy)  is  the  Mueller  matrix  for  the  scatterer  as  obtained  from 
the  coupled-dipole  method.  To  determine  the  amount  of  scattered  light  that  is 
horizontally  polarized,  the  Mueller  matrix  for  a  horizontal  polarizer  is  required 
once  more.  The  product  of  these  matrices  gives  the  Stokes  vector  for  scattered 
light.  The  I,  component  of  the  resulting  Stokes  vector  represents  the  irradiance 
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of  horizontally  polarized  light  scattered  by  a  particle  that  was  illuminated  by 
horizontally  polarized  light.  This  irradiance  will  be  referred  to  as  SHhi  and  its 
analog  for  vertically  polarized  light  will  be  called  Sw.  The  product  of  the 
appropriate  Mueller  matrices  yields  the  following  equations: 


$HH  "  ^ll+^12+^2l+^22^ 

Syy  ”  ~ O^i i — ^  12~*^21  +l^22^ 


(4.9) 


Bickel  and  Bailey  (1985)  provide  an  illuminating  characterization  of  scattered 
light  in  the  context  of  Stokes  vectors  and  Mueller  matrices. 


Cross  Sections 


If  a  particle  is  interposed  between  a  source  of  electromagnetic  radiation 
and  a  detector,  the  power  received  at  the  detector  decreases.  The  decrease  in 
power  or  extinction  of  the  incident  wave  is  due  to  scattering  and  absorption  by 
the  particle  and  depends  on  the  particle’s  shape,  size,  and  relative  complex 
refractive  index  as  well  as  the  frequency  and  polarization  state  of  the  incident 
plane  wave.  Characterization  of  the  extinction  is  usually  denoted  in  terms  of 
particle  extinction  C„„  scattering  C^,  and  absorption  Cabs  cross  sections.  For 
plane  waves,  each  cross  section  can  be  determined  independently,  and 
agreement  can  be  tested  by  considering  energy  conservation 
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c  -  c  .  +  C  .  (4.10) 

ext  abs  sea  '  v  7 

The  calculation  of  these  cross  sections  in  the  coupled-dipole  method  is  now 
briefly  described. 


Extinction 


The  extinction  cross  section  is  calculated  using  the  optical  theorem 
which  indicates  that  extinction  depends  only  on  the  scattering  amplitude  in  the 
forward  direction.  This  may  seem  counterintuitive  since  extinction  accounts 
for  both  absorption  in  the  particle  and  scattering  in  all  directions  by  the 
particle.  For  plane  waves,  the  extinction  cross  section  is  determined  by 
calculating  the  work  done  by  the  incident  electric  field 


ext 


(4.11) 


where  pj  is  the  dipole  moment  and  Ej  is  the  incident  electric  field  at  the  jth 
dipole  as  calculated  by  the  coupled-dipole  method,  *  indicates  the  complex 
conjugate,  and  Im  signifies  that  only  the  imaginary  part  of  the  argument  be 


evaluated. 
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Scattering 


The  scattering  cross  section  is  determined  by  first  calculating  the  time- 
averaged  power  radiated  per  unit  solid  angle  (Jackson,  1975) 


- -  E  l(«x^)x«|2. 


(4.12) 


where  dC^/dQ  is  sometimes  called  the  differential  scattering  cross  section. 
This  quantity  represents  the  amount  of  light  scattered  into  a  unit  solid  angle  in 
direction  n.  Simpson’s  rule  is  used  to  integrate  over  the  total  solid  angle  to 
yield  the  scattering  cross  section.  For  a  good  estimate  of  scattering  cross 
section,  it  is  necessary  to  sum  over  a  minimum  of  33  theta  ( e )  and  12  phi  (0) 
values,  where  dQ  =  sin0dtfd0  (Draine,  1988). 

The  average  cosine  of  the  scattering  angle,  or  the  asymmetry  parameter 
g,  can  be  obtained  while  calculating  the  scattering  cross  section  from 
equation  4.12  as 


g  =  (cosQ)  - - - —  f dQn-k  \'Z[p.-n(ri'p)]e  (4.13) 

C  |E.|2J  M  1  1 

sca\  1 1 

The  asymmetry  parameter  ranges  in  value  from  -1  to  +1.  For  a  particle  whose 
scattered  power  density  is  symmetric  about  the  scattering  angle  of  90°,  g 
vanishes  identically.  If  the  forward  scattering  is  larger  (smaller)  than 
backscattering,  then  the  asymmetry  parameter  is  positive  (negative). 
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Absorption 


The  absorption  cross  section  is  determined  by  summing  the  power 
absorbed  by  all  the  dipoles 


4izk0  JL  pX  p 

- l Imp:'  —  -  — k3p . 'p: 

|jEi|-  ;■-!  L  V a  /  J  * 


(4.14) 


The  first  term  on  the  right  is  the  expression  used  by  Purcell  and  Pennypacker 
(1973).  Draine  (1988)  added  the  second  term  (radiative  reaction  term)  to 
account  for  effect  of  the  radiation  on  the  motion  of  the  dipoles  as  prescribed 
by  Jackson  (1975,  Chap  17).  Including  the  second  term  improves  calculation  of 
Ca5s  when  comparing  absorption  by  spheres  with  Mie  theory. 

The  accuracy  of  these  equations  that  determine  the  three  cross  sections 
can  be  shown  by  comparing  C +  Cabs  with  Cext.  Table  4.2  shows  these  values 
for  several  spheres  with  different  refractive  indices. 


Buckscattering 


The  backscattering  cross  section  is  calculated  using  the  traditional 


definition  of  the  radar  backscattering  cross  section  (see  Bohren  and  Huffman, 
1983,  Sect  4.6) 
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Table  4.2  Results  from  the  three  independent  equations  used  to 
determine  CKi,  C^,  and  Cext.  Conservation  of  energy 
dictates  that  Cxa  +  Cabl  =  Cext.  Values  in  table  show 
excellent  agreement.  Values  calculated  by  Mie  theory  are 
also  presented:  Cext(Mie).  Size  parameter  of  the  modeled 
sphere  is  1.6.  Cross  sections  have  units  of  area. 


Refractive  index 

c 

^sca 

Cabs 

Cabs  +  Oca 

cext 

Ce«(Mie) 

1.14  +  i0.26 

5.87 

27.00 

32.87 

32.88 

33.23 

1.33  +  i0.05 

11.09 

7.65 

18.74 

18.75 

19.22 

1.55  +  iO.OOS 

33.84 

1.00 

34.84 

34.84 

37.32 

1.39  +  i0.42 

17.77 

40.14 

57.91 

57.91 

57.99 

1.7  +  iO.l 

47.78 

19.23 

67.01 

67.02 

72.43 

1.9  +  i0.0004 

108.78 

0.13 

108.91 

108.92 

121.40 

2.5  +  il.4 

50.03 

61.49 

111.52 

111.53 

97.85 

3.5  +  i2.05 

52.37 

58.29 

110.66 

110.67 

94.28 

3.0  +  i4.0 

53.59 

56.11 

109.70 

109.70 

92.97 
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Ctel  -  4~(l52(180°)l2+  |5,(180°)P),  (4.15) 

At 

where  S2(180°)  and  St(180°)  are  elements  from  the  amplitude  scattering 
matrix. 


The  cross  section  efficiencies  Q  are  calculated  by  dividing  the  respective 


cross  section  by  the  geometric  cross  section  of  the  particle  G  normal  to  the 
incident  wave.  For  example,  the  extinction  efficiency  is  determined  by 
Q„  =  C„/G. 


Chapter  5 

BACKSCATTERING  CALCULATIONS 
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Light  scattering  is  important  in  many  fields  of  science:  for  example, 
inferring  properties  of  interstellar  dust  (Draine,  1988),  locating  dust  particles 
on  computer  chips  (Taubenblatt,  1990),  studying  molecules  by  analyzing  the 
polarization  state  of  scattered  light  (Harris  and  McClain,  1985).  In  remote 
sensing  of  the  atmosphere,  often  the  only  information  available  is  the 
backscattered  signal.  The  objective  of  this  chapter  is  to  investigate  the 
advantages  and  disadvantages  as  well  as  the  limitations  of  using  the  coupled- 
dipole  method  for  determining  backscattering  by  hydrometeors. 

Although  most  solid  particles  in  nature  are  not  spherical,  it  has  become 
too  often  the  practice  to  model  these  particles  as  spheres  using  Mie  theory. 
Thus,  the  modeler  avoids  using  a  more  computer  intensive,  but  in  many  cases 
more  appropriate,  model.  Unfortunately,  the  assumption  that  randomly 
oriented,  nonspherical  particles  scatter  like  equal-volume  spheres  has  become 
the  mindset  of  many  modelers.  Methods  for  calculating  scattering  by  arbitrary 
particles  are  often  used  to  identify  circumstances  under  which  the  practice  of 
using  equal-volume  spheres  is  least  likely  to  lead  to  unacceptable  errors,  but  it 
is  done  in  a  fashion  that  appears  to  condone  this  practice  rather  than  condemn 
it.  In  this  chapter,  calculations  for  arbitrary  particles  will  be  compared  with 
those  for  equivalent-volume  spheres. 
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The  first  topic  in  this  chapter  is  the  shape  sensitivity  of  the 
backscattered  signal.  Next,  scattering  results  from  modeling  spherical  particles 
using  the  coupled-dipole  method  are  compared  with  Mie  theory,  and  emphasis 
is  placed  on  how  accuracy  depends  on  size  parameter.  This  is  followed  by  a 
brief  discussion  of  how  well  a  collection  of  dipoles  on  a  cubic  lattice  can 
represent  a  solid  sphere,  and  whether  small  variations  of  shape  affect 
backscattering  calculations.  Next,  backscattering  by  equal- volume  particles  of 
different  shapes  are  compared.  With  information  from  the  preceding  analyses, 
backscattering  of  94  GHz  radar  by  ice  spheres  and  crystals  is  investigated. 
First,  scattering  results  from  modeling  spheres  with  the  coupled-dipole  method 
are  compared  with  Mie  theory  to  determine  the  accuracy  of  the  model  for  the 
refractive  index  of  ice  at  94  GHz:  1.878  +  i4.76x  10-4.  Finally,  hexagonal 
plates  and  columns  are  modeled.  The  impetus  is  to  identify  the  advantages  of 
using  linearly  polarized  radar  to  distinguish  shape  and  size  dependencies  in 
backscattering  by  ice  crystals. 

Shape  and  Size  Dependency  of  Backscattering 

Bad  data  can  be  made  to  look  good  with  a  suitable  choice  of  statistics. 
In  a  similar  fashion,  numerical  models  can  appear  credible  if  only  forward 
scattering  calculations  are  compared  with  analytic  solutions.  The  coupled- 
dipole  method  can  be  used  to  demonstrate  why  scattering  is  more  critical  in 
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the  backward  direction.  The  scattering  light  is  the  sum  of  all  the  radiated 
waves  from  the  dipolar  subunits,  the  intensity  of  the  scattered  light  depends  on 
the  phase  relation  among  all  the  radiated  waves.  The  following  argument 
appears  in  Bohren  and  Singham  (1990). 

The  scattered  field  is  the  sum  of  the  radiated  waves  by  all  the  dipoles 
that  represent  the  particle.  Consider  the  phase  difference  a<£  of  two  waves 
scattered  by  two  separate  dipoles  excited  by  the  same  wave  as  shown  in 
Figure  5.1.  Interaction  between  the  dipoles  is  ignored.  The  phase  difference 
depends  on  the  distance  between  the  dipoles  d  (in  the  direction  of  wave 
propagation),  the  wavelength  of  incident  light,  and  the  angle  between  the 
incident  and  scattered  waves  6: 


«» 
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For  large  particles  when  the  phase  difference  of  the  incident  light  across 
the  particle  cannot  be  neglected,  scattering  is  affected  by  interference  among 
the  radiated  waves,  which  can  be  accounted  for  by  the  dipole  interactions. 

Since  the  coupled-dipole  method  accounts  for  the  interactions,  it  is  expected  to 
accurately  compute  backscattering.  This  is  illustrated  in  Figure  5.2  where 
scattering  results  using  the  coupled-dipole  method  are  compared  with  Mie 
theory.  The  sphere  modeled  has  a  refractive  index  of  1.33  +  i0.05  and  a  si2e 
parameter  of  1.60.  The  array  used  in  the  coupled-dipole  method  contains  461 
dipoles;  computations  were  made  using  the  scattering-order  technique.  The 
solid  line  represents  Mie  calculations.  The  dotted  line  is  from  the  coupled- 
dipole  method  after  the  first  iteration,  which  is  equivalent  to  ignoring  the 
dipole  interactions;  the  dashed  line  represents  a  fully  converged  solution 
(14  iterations).  Backscattering  calculations  are  more  sensitive  to  dipole 
interactions  than  forward  scattering  calculations.  Thus,  when  judging  model 
capability  it  is  better  to  compare  backscattering  calculations,  since  forward 
scattering  is  least  sensitive  to  shortcomings  in  the  model. 

The  difficulty  in  accurately  calculating  the  backscattered  signal  increases 
with  size  parameter.  Figure  5.3  shows  forward  scattering  (dashed  line)  and 
backscattering  (solid  line)  as  calculated  by  Mie  theory  for  a  sphere  with 
refractive  index  of  1.33  +  i0.05  and  illuminated  by  unpolarized  light.  For  size 
parameters  less  than  about  0.5  both  scattered  signals  vary  with  the  sixth  power 
of  particle  radius,  which  is  indicative  of  particles  in  the  Rayleigh  regime 
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Figure  5.2  Importance  of  dipole  interactions  for  backscattering  calculations. 

A  sphere  with  refractive  index  133  +  i0.05  and  size  parameter  of 
1.60  was  modeled.  Solid  line  is  from  Mie  theory.  The  other 
lines  are  from  the  coupled-dipole  method  using  the  scattering- 
order  technique.  The  dotted  line  represents  scattering  values 
with  no  dipole  interactions;  the  dashed  line  represents  a  fully 
converged  solution  that  includes  dipole  interactions. 
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Figure  53  Forward  scattering  and  backscattering  by  spheres  of  increasing 
size  parameter.  Spheres  with  refractive  index  133  +  i0.05  and 
size  parameters  ranging  from  0.1  to  50  were  modeled.  The 
dashed  line  represents  forward  scattering;  the  solid  line 
represents  backscattering. 
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(negligible  phase  shift  of  the  incident  wave  across  the  particle).  For  larger 
particles,  the  effects  of  interference  cause  backscattering  to  fluctuate  wildly; 
hence,  small  errors  in  particle  size  further  reduce  computational  accuracy.  The 
values  in  Figure  5.3  are  obtained  from  Mie  theory,  which  is  a  solution  for  a 
sphere  with  a  well-defined  radius.  In  the  coupled-dipole  method  the 
representation  of  arbitrary  particles  by  arrays  of  dipoles  can  make  computing 
backscattering  less  precise.  Limitations  are  now  examined. 


The  ability  of  the  coupled-dipole  method  to  calculate  backscattering 
depends  on  the  particle’s  size  and  refractive  index  and  how  well  the  array  of 
dipoles  represents  the  particle’s  shape.  Dependence  on  refractive  index  is 
addressed  in  Appendix  A.  Table  A.l  shows  that  accuracy  of  scattering 
calculations  from  the  coupled-dipole  method  decreases  with  increasing 
refractive  index.  This  is  due  to  the  influence  of  higher  order  multipoles.  In 
this  section  the  effect  of  small  variations  in  particle  shape  and  of  increasing  the 
particle  size  will  be  examined. 

Changing  the  particle  size  affects  the  accuracy  of  the  scattering 
calculations  using  the  coupled-dipole  method.  When  the  particle  is  small 
compared  with  the  wavelength  (size  parameter  less  than  0.5)  the  phase 
difference  of  the  incident  light  across  the  particle  is  negligible.  In  turn,  the 
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dipoles  are  also  small  compared  with  the  wavelength.  Increasing  the  particle 
size  requires  either  increasing  the  radius  or  number  of  dipoles.  The  latter 
option  is  permissible  as  long  as  the  radius  of  the  dipoles  satisfies  the  criterion 
that  koa  |  m  |  <  0.5,  and  this  criterion  is  not  violated  in  the  following  analysis. 

A  sphere  with  a  refractive  index  of  1.33  +  i0.05  was  modeled;  its  size 
parameter  was  varied  from  1.0  to  2.0  by  increasing  the  radius  of  the  dipoles. 
The  sphere  was  represented  by  three  separate  arrays:  251,  304,  and  461 
dipoles.  The  results  in  Figure  5.4  show  good  agreement  with  Mie  theory  until 
the  size  parameter  exceeds  1.6.  It  is  at  this  size  parameter  that  the 
backscattered  signal  reaches  its  first  minimum  value  (as  shown  in  Figure  5.3). 
This  poor  agreement  with  Mie  theory  at  size  parameters  above  that  which 
corresponds  to  the  backscattering  minimum  also  occurs  for  other  refractive 
indices.  These  results  are  the  same  whether  using  the  matrix  inversion  or 
scattering-order  technique.  To  a  lesser  degree,  but  still  occurring  at  the 
backscattering  minimum,  forward  scattering  calculations  lose  accuracy  as  shown 
in  Figure  A.2.  Thus,  for  size  parameters  larger  than  that  which  corresponds  to 
the  first  backscattering  minimum,  scattering  calculations  from  the  coupled- 
dipole  method  lose  accuracy. 

This  loss  of  accuracy  can  be  overcome  somewhat  by  increasing  the 
number  of  dipoles  in  the  array  while  decreasing  their  size.  Figure  5.5  shows 
the  backscattering  for  a  sphere  with  refractive  index  1.33  +  iO.l  represented  by 
arrays  of  461,  739,  1064,  and  2320  dipoles;  results  from  Mie  theory  are 
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Figure  5.4  Comparison  of  coupled-dipole  method  with  Mie  theory  for 

spheres  of  increasing  size  parameter.  A  sphere  with  refractive 
index  133  +  i0.05  was  modeled.  The  particle  size  was  increased 
by  increasing  the  size  of  the  dipoles.  The  solid  line  represents 
scattering  by  a  461  dipolar  array,  the  dotted  line  by  a  304  dipolar 
array,  and  the  dashed  line  by  a  251  dipolar  array. 
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Figure  5.5 


Improvement  of  backscattering  by  increasing  the  number  of 
dipoles.  A  sphere  with  refractive  index  133  +  iO.l  and  a  size 
parameter  1.60  was  model.  The  solid  line  is  Mie  theory.  The 
short-dashed  line  (-  -  -)  represents  scattering  by  a  461  dipolar 

array,  the  multiple-dashed  line  ( - )  by  a  739  dipolar  array, 

the  dotted  line  by  a  1064  dipolar  array,  and  the  long-dashed  line 
( - )  by  a  2320  dipolar  array. 
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depicted  by  the  solid  line.  In  general,  backscattering  calculations  improve  by 
increasing  the  number  of  dipoles;  computation  time  always  increases. 

Increasing  the  number  of  dipoles  in  the  array  produces  subtle  changes  in 
the  shape  of  the  array.  The  question  of  how  well  an  array  of  dipoles  on  a 
simple  cubic  lattice  represents  a  spherical  particle  is  addressed  in  Appendix  A 
beginning  on  page  116.  Draine  (1988)  presented  a  formula  for  estimating  the 
sphericity  of  an  array  of  dipoles  that  represents  a  sphere.  Table  A.2  shows 
that  small  variations  in  sphericity  may  affect  scattering  when  shape  dependency 
is  critical.  As  discussed  in  the  previous  section,  backscattering  is  very  shape 
dependent;  therefore,  if  these  small  variations  of  sphericity  are  important  they 
should  affect  backscattering  calculations. 

To  check  how  backscattering  depends  on  shape  variations  associated 
with  representing  a  sphere  by  a  lattice  of  dipoles,  two  spheres  with  refractive 
index  1.33  +  iO.l  were  modeled;  one  had  a  size  parameter  of  0.84  and  the 
other  1.60.  Small  variations  in  sphericity  were  obtained  by  representing  each 
sphere  with  13  arrays  that  varied  in  size  from  280  to  515  dipoles.  The 
effective  radius  of  each  array  was  held  constant  by  varying  the  radius  of  the 
dipoles.  Thus,  the  main  difference  between  the  arrays  was  the  specific  shape 
of  each  array.  By  Draine’s  criterion  (1988),  a  perfect  sphere  has  a  radius  of 
gyration  of  1,  and  the  value  for  radius  of  gyration  increases  at  the  particle 
departs  from  sphericity.  Backscattering  results  are  shown  in  Table  5.1.  For 
the  smaller  sphere  the  backscattering  values  all  agree  with  Mie  theory  to 
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Table  5.1  Effect  of  small  variations  in  shape  on  backscattering  calculations. 

Two  spheres  with  refractive  index  of  1.33  +  iO.l  were  modeled; 
their  size  parameters  are  0.84  and  1.60.  N  is  the  number  of 
dipoles  in  the  arrays,  ROG  is  the  radius  of  gyration  of  the  array 
(Draine,  1988).  ROG  is  an  estimate  of  the  sphericity  of  the 
array--an  exact  sphere  has  a  value  of  1,  deviation  from  1 
indicates  nonsphericity.  Rows  are  arranged  by  increasing  ROG. 
The  percent  difference  calculation  is: 

[Stl(180°)/S11(180°;  Mie)]-1.  For  the  smaller  sphere,  effects  of 
small  variations  in  shape  are  not  critical  as  for  the  larger  sphere. 


N 

ROG 

%  diff  when 
x  =  0.84 

%  diff  when 
x  =  1.60 

461 

1.0006 

0.6% 

-2.6% 

304 

1.0012 

0.6 

9.1 

280 

1.0018 

1.3 

15.9 

437 

1.0019 

1.1 

17.0 

485 

1.0020 

0.1 

-10.3 

480 

1.0027 

1.6 

28.8 

432 

1.0031 

0.6 

-2.1 

515 

1.0034 

0.6 

-2.4 

360 

1.0035 

0.2 

8.2 

389 

1.0035 

1.4 

18.0 

365 

1.0039 

0.9 

8.1 

341 

1.0048 

1.5 

25.1 

305 

1.0060 

1.1 

8.2 

66 


within  2%;  for  the  larger  sphere,  less  than  one-quarter  of  the  values  agree  to 
within  3%  while  nearly  half  disagree  by  more  than  10%.  The  poorer 
agreement  for  the  larger  particle  is  because  its  size  is  beyond  the  first 
backscattering  minimum  as  discussed  on  page  58.  At  this  size  parameter 
backscattering  becomes  extremely  shape  dependent.  A  relation  between 
increased  sphericity  and  improved  backscattering  calculations  is  not  observed 
in  Table  5.1.  Thus,  Draine’s  criterion  for  radius  of  gyration  does  not  generally 
aid  in  choosing  an  array  that  best  represents  a  sphere.  This  criterion  may  have 
merit  for  calculations  in  the  forward  direction;  however,  for  the  more  shape- 
dependent  backward  direction  it  provides  no  guidance  for  selecting  arrays. 

For  particles  that  are  larger  than  that  which  produces  the  first 
backscattering  minimum,  scattering  calculations  using  the  coupled-dipole 
method  become  suspect.  Increasing  the  number  of  dipoles  in  the  array 
generally  improves  the  accuracy  of  these  calculations. 

Backscattering  by  Equal-Volume  Particles 

It  was  shown  in  the  earlier  sections  of  this  chapter  that  accuracy  of 
backscattering  calculations  using  the  coupled-dipole  method  is  dependent  on 
particle  size  and  shape.  This  last  sensitivity  study  examines  whether  the 
overall  particle  shape  strongly  affects  backscattering  calculations.  The  equal- 
volume  particles  modeled  here  are  a  prolate  spheroid,  rectangular  solid, 
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cylinder,  and  hexagonal  column;  all  have  an  aspect  ratio  (ratio  of  length  to 
diameter)  of  5:1  and  size  parameter  of  1.0.  An  equal-volume  sphere  is 
included.  Three  refractive  indices  chosen  are:  1.5  +  iO.OO,  1.5  +  i0.05,  and 
1.5  +  i0.50;  these  values  are  representative  of  ice  crystals  in  the  wavelength 
range  from  6  to  15  /urn  (peak  outgoing  terrestrial  radiation)  (Rusk  et  al.,  1971). 

Scattering  calculations  were  made  with  these  particles  in  orientations 
that  simulate  natural  conditions.  Except  for  very  large  sizes  and  under 
turbulent  conditions,  ice  crystals  are  expected  to  fall  without  tumbling  (Cho  et 
al.,  1981;  Pruppacher  and  Klett,  1980);  therefore,  the  long  axis  of  the  particles 
modeled  here  remain  horizontal  (in  the  x-y  plane).  The  particles  were  allowed 
to  rotate  in  the  horizontal  plane,  but  not  spin  about  their  longitudinal  axis.  To 
account  for  not  spinning  about  the  longitudinal  axis,  the  particle  was  initially 
oriented  in  the  position  that  yielded  the  average  backscattering  value.  For 
instance,  the  rectangular  solid  was  rotated  45°  and  the  hexagonal  crystal  by 
22.5°;  an  initial  rotation  did  not  significantly  alter  backscattering  by  the  prolate 
spheroid  or  the  cylinder  because  of  their  symmetry.  Several  states  of 
polarization  of  the  incident  wave  were  analyzed:  Sn,  Shh,  and  Sw  (see  page  47 
for  a  description  of  these).  Modeling  results  are  shown  in  Figures  5.6  through 
5.9. 

Figure  5.6  shows  Stl(180°)  for  zenith  angles  ranging  from  0°  to  90°  (the 
angle  is  defined  with  the  zenith  at  0°  and  the  horizon  at  90°)  for  the  five 
particles  with  refractive  index  1.5  +  iO.OO.  The  horizontal  line  represents 


68 


Figure  5.6  Backs  cattering  of  unpolarized  light  by  equal-volume  particles 

with  aspect  ratios  of  5:1  and  refractive  index  13  +  iO.OO.  The 
solid  horizontal  line  represents  scattering  by  an  equal-volume 
sphere.  The  other  solid  line  represents  scattering  by  a  hexagonal 
column,  the  short-dashed  line  by  a  prolate  spheroid,  the  dotted 
line  by  a  cylinder,  the  multiple-dashed  line  by  a  rectangular  solid. 
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Figure  5.7  Backscattering  of  horizontally  polarized  light  by  equal-volume 
particles  with  aspect  ratios  of  5:1  and  refractive  index 
1.5  +  iO.OO.  The  solid  horizontal  line  represents  scattering  by  an 
equal-volume  sphere.  The  other  solid  line  represents  scattering 
by  a  hexagonal  column,  the  short-dashed  line  by  a  prolate 
spheroid,  the  dotted  line  by  a  cylinder,  the  multiple-dashed  line 
by  a  rectangular  solid. 
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Figure  5.8  Backscattering  of  vertically  polarized  light  by  equal-volume 
particles  with  aspect  ratios  of  5:1  and  refractive  index 
1.5  +  iO.OO.  The  solid  horizontal  line  represents  scattering  by  an 
equal-volume  sphere.  The  other  solid  line  represents  scattering 
by  a  hexagonal  column,  the  short-dashed  line  by  a  prolate 
spheroid,  the  dotted  line  by  a  cylinder,  the  multiple-dashed  line 
by  a  rectangular  solid. 
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Figure  5.9  Backscattering  of  unpolarized  light  by  equal-volume  particles 

with  aspect  ratios  of  5:1  and  refractive  index  of  1.5  +  iO.50.  The 
solid  horizontal  line  represents  scattering  by  an  equal-volume 
sphere  The  other  solid  line  represents  scattering  by  a  hexagonal 
column,  the  short-dashed  line  by  a  prolate  spheroid,  the  dotted 
line  by  a  cylinder,  the  multiple-dashed  line  by  a  rectangular  solid. 
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backscattering  by  the  equal-volume  sphere  (no  angle  dependence).  The 
remaining  four  particles  display  a  similar  relationship  between  zenith  angle  and 
backscattering.  These  results  indicate  that  the  overall  particle  shape  is  more 
important  that  subtle  differences  in  the  specific  cross-section  shape.  Thus,  for 
particles  with  aspect  ratios  of  at  least  5:1,  analytic  or  numerical  methods  which 
treat  elongated  spheroids  may  be  used  to  calculate  backscattering  by  ice 
columns  with  only  minor  compromise  to  model  accuracy  (e.g.  Yeh  et  al.,  1982). 

Figures  5.7  and  5.8  show  8^(180°)  and  Sw(180°)  for  the  same  particles 
and  zenith  angles  as  in  Figure  5.6.  At  zenith  angle  of  0°,  corresponding  values 
of  Shh(180°)  and  Sw(180°)  are  equal  because  at  this  angle  the  particles  appear 
symmetric.  With  increasing  zenith  angle  (looking  down  towards  the  horizon), 
the  values  of  SHH(180°)  decrease  slower  than  those  for  Sw(180°).  This  too  is 
expected  because  the  horizontally  polarized  light  excites  the  longer  dimension 
of  the  particle.  This  difference  between  Shh(180°)  and  Sw(180°)  is  important, 
and  it  will  be  discussed  later  as  a  means  to  infer  the  size  distribution  of  ice 
crystals  by  radar. 

When  the  refractive  index  of  these  particles  is  changed  to  1.5  +  i0.05, 
the  backscattering  calculations  are  similar  to  those  shown  in  Figures  5.6 
through  5.8  and  are  not  presented.  Sn(180°)  values  for  these  particles  with 
refractive  index  1.5  +  i0.50  are  shown  in  Figure  5.9.  Backscattering  for  all 
three  refractive  indices  displays  similar  angle  dependencies.  The  only 
differences  of  note  are  the  higher  values  when  the  refractive  index  is 
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differences  of  note  are  the  higher  values  when  the  refractive  index  is 
1.5  +  i0.50.  This  is  a  consequence  of  the  smaller  penetration  depth  due  to  the 
larger  imaginary  part  of  the  refractive  index.  Values  for  Shh(180°)  and 
Sw(180  )  at  this  refractive  index  behave  similarly  to  those  in  Figures  5.7 
and  5.8. 

The  mean  value  of  Sn(180°)  (averaged  over  all  zenith  angles)  for  the 
four  elongated  particles  is  similar  to  the  backscattering  by  the  equal-volume 
sphere.  The  percent  difference  between  the  mean  values  of  Su(180°)  for  the 
elongated  particles  and  the  backscattering  value  for  the  equal-volume  sphere 
increases  with  increasing  imaginary  part  of  the  refractive  index  (i.e.,  3.6%, 

3.8%,  and  9.1%  for  these  three  refractive  indices).  Mugnai  and  Wiscombe 
(1980)  noted  that  absorption  always  improves  the  agreement  of  backscattering 
by  equivalent  spheres;  these  calculations  indicate  the  opposite.  However,  they 
modeled  particles  that  had  only  a  slight  deviation  (10%)  from  sphericity. 

The  backscattering  calculations  reported  in  this  chapter  show  size  and 
shape  dependence.  In  addition,  it  is  apparent  that  backscattering  calculations 
using  the  coupled-dipole  method  lose  accuracy  when  the  particle  size  is  larger 
than  that  associated  with  the  first  backscattering  minimum.  It  also  appears  that 
for  elongated  panicles,  small  variations  in  the  particle’s  cross  section  are 
inconsequential.  The  last  section  addresses  the  size  and  shape  dependence  of 
the  first  backscattering  minimum  of  ice  crystals  at  94  GHz. 
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Backscattering  of  94  GHz  Radar  bv  Ice  Crystals 

The  coupled-dipole  method  is  applicable  to  studying  backscattering  of 
94  GHz  radar  by  ice  crystals  for  several  reasons.  First,  the  coupled-dipole 
method  can  calculate  scattering  by  ice  crystals.  Second,  the  wavelength  of  this 
radar  is  3 2  mm;  hence,  the  size  parameter  (or  approximately  the  number  of 
dipoles  N  needed  to  represent  the  particle)  is  not  prohibitively  large  when 
modeling  particles  with  linear  dimensions  of  the  order  of  one  to  five 
millimeters.  Third,  with  N  <  300  the  matrix  inversion  technique  can  be  used 
to  rotate  the  particle  through  many  orientations.  Lastly,  the  refractive  index  of 
ice  at  this  frequency  is  not  too  large  to  prevent  use  of  the  scattering-order 
technique  for  calculating  scattering  by  larger  ice  crystals  (when  N  >  300). 

Of  the  models  that  calculate  scattering  by  arbitrary  particles,  the 
coupled-dipole  method  is  best  suited  for  studyting  ice  crystals.  For  example, 
the  T-matrix  method  calculates  scattering  by  rotationally  symmetric  particle 
with  low  aspect  ratios;  the  ice  crystals  modeled  here  do  not  fit  that  description. 

The  relationship  between  particle  shape  and  the  first  backscattering 
minimum  is  one  of  the  interests  for  studying  backscattering  at  94  GHz.  For 
falling  raindrops,  a  backscattering  minimum  occurring  at  radius  0.83  mm 
provides  an  opportunity  to  remotely  identify  spherical  raindrop  size  (Lhermitte, 
1988).  This  identification  is  possible  because  the  terminal  velocity  of  raindrops 
increases  with  drop  radius  while  backscattering  by  raindrops  is  at  a  local 
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minimum  for  drop  radius  0.83  mm.  Doppler  radars  measure  the  line-of-sight 
velocity  of  its  targets  (raindrops).  When  the  radar  is  looking  above  the  horizon 
(elevation  angle  greater  than  0°),  the  Doppler  shift  is  a  result  of  both  the 
horizontal  and  vertical  forces  acting  on  the  target.  The  Doppler  shift  of  a 
falling  raindrop  increases  with  drop  size  because  terminal  velocity  varies  with 
drop  size.  If  the  Marshall-Palmer  (1948)  raindrop  size  distribution  is  assumed, 
the  relationship  between  the  backscattering  for  each  raindrop  size  and  its 
respective  Doppler  shift  can  be  determined.  The  crucial  point  is  that  the 
backscattering  minima  are  also  present  in  the  Doppler  shift  spectrum.  Since 
the  terminal  velocity  of  raindrops  with  radius  0.83  mm  is  known,  the  vertical 
component  of  the  Doppler  shift  at  that  minimum  can  be  determined.  By 
removing  the  vertical  component  of  the  Doppler  shift,  the  horizontal 
component  remains.  From  the  horizontal  component  one  obtains  a  better 
estimate  of  the  clear  air  velocity.  This  information  about  air  motion  may  help 
to  better  understand  cloud  dynamics,  cloud  structure,  and  precipitation  process. 

This  procedure  is  based  on  identifying  the  backscattering  minimum 
which  occurs  for  spheres.  If  a  backscattering  minimum  is  present  for  a  size 
distribution  of  ice  crystals  this  procedure  could  increase  our  understanding  of 
ice  clouds.  Since  ice  crystals  are  nonspherical,  Mie  theory  is  not  an  applicable 
modeling  tool;  therefore,  an  alternative  model  such  as  the  coupled-dipole 
method  is  required.  Before  calculating  the  backscattering  of  94  GHz  radar  by 
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ice  crystals,  the  accuracy  of  the  coupled-dipole  method  will  be  estimated  by 
first  comparing  scattering  results  for  ice  spheres  with  Mie  theory. 

Ice  Spheres 

The  scattering-order  technique  was  used  to  calculate  backscattering  by 
ice  spheres,  and  the  results  were  compared  with  Mie  theory.  The  refractive 
index  of  ice  (0°  C)  at  94  GHz  is  1.878  +  i4.76  x  10-4.  The  number  of  dipoles 
in  the  arrays  used  to  represent  the  sphere  ranged  from  461  for  the  smallest 
sphere  to  2320  for  the  largest  sphere.  The  number  of  dipoles  was  increased  to 
keep  the  dipole  radius  small.  This  improves  the  accuracy  of  the  calculations 
above  the  backscattering  minimum,  which  occurs  at  a  radius  of  0.73  mm  for 
ice.  The  results  shown  in  Figure  5.10  indicate  that  calculations  by  the  coupled- 
dipole  method  do  fairly  well  at  identifying  the  correct  size  parameter  and 
magnitude  of  the  backscattering  response  at  the  first  minimum.  The  final  step 
is  to  calculate  backscattering  by  ice  columns  and  plates. 

Hexagonal  Columns  and  Plates 

The  two  ice  crystals  selected  to  be  modeled  are  a  hexagonal  column  and 
plate  with  aspect  ratios  3.5:1  and  3.2:1,  respectively.  Hexagonal  prismatic  is 
the  basic  shape  of  ice  crystals,  although  laboratory  observations  have  revealed 
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a  few  other  shapes  (Pmppacher  and  Klett,  1980,  Chap  2).  The  aspect  ratio  of 
5:1  for  the  column  would  have  been  a  better  representation  based  on 
observations,  but  computer  memory  requirements  restricted  the  size. 

The  matrix  inversion  technique  was  used  to  calculate  backscattering  by 
the  ice  crystals.  The  first  set  of  calculations  were  made  to  dispel  the  notion 
that  randomly  oriented  particles  always  scatter  like  equal-volume  spheres. 

Backscattering  by  the  ice  column  and  plate  was  calculated  with  the 
particles  positioned  in  100  orientations.  The  values  for  the  column,  plate,  and 
sphere  are  shown  in  Figure  5.11.  For  size  parameters  less  than  0.8  the 
difference  between  the  backscattered  signals  is  not  too  great.  (This  was 
observed  in  Figures  5.6  through  5.9  where  the  average  backscattered  signal 
from  the  randomly  oriented  particles  with  a  size  parameter  of  1.0  was 
approximately  equal  to  the  backscatter  from  an  equal-volume  sphere.)  At  size 
parameters  larger  than  1.0  the  backscattering  signal  of  the  sphere  does  not 
coincide  with  the  minima  for  the  other  particles.  Beyond  this  particle  size,  the 
backscattered  signals  differ  by  as  much  as  a  factor  of  six.  Thus,  at  small  size 
parameters  and  for  particles  of  refractive  indices  close  to  unity,  using  an  equal- 
volume  sphere  in  place  of  randomly  oriented  particles  may  not  be  a  bad 
assumption  when  calculating  orientationally  averaged  backscattered  signals; 
however,  at  larger  size  parameters  the  backscattering  becomes  very  shape 
dependent,  randomly  oriented  particles  might  no  longer  scatter  like  spheres. 
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Figure  5.11  Backscattering  of  unpolarized  light  by  randomly-oriented,  equal- 
volume  particles  of  increasing  size  parameter  and  with  refractive 
index  1.878  +  iO.OGOS.  Calculations  were  made  with  the  coupled- 
dipole  method.  The  solid  line  represents  a  sphere,  the  dashed- 
line  represents  a  hexagonal  column  with  aspect  ratio  of  3-5:1,  and 
the  dotted  line  represents  a  hexagonal  plate  with  aspect  ratio  of 
32:1. 
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The  next  set  of  calculations  was  made  to  determine  a  shape  and  size 
dependence  of  the  first  backscatter  minimum  for  ice.  For  subsequent  modeling 
the  orientation  of  the  crystals  was  assumed  to  be  identical  to  that  in  the 
preceding  section  on  backscattering  by  equal-volume  particles:  the  particles  do 
not  tumble,  and  they  are  allowed  to  rotate  within  the  horizontal  plane. 

Figures  5.12  through  5.16  show  the  dependence  of  backscattering  by  the 
hexagonal  column  (aspect  ratio  of  3.5:1)  on  size  parameter,  which  is  based  on 
that  of  an  equal-volume  sphere.  The  values  for  size  parameter  range  from  1.0 
to  2.2  and  correspond  to  column  lengths  of  2.1  to  4.6  mm.  The  five  lines  in 
the  figures  represent  zenith  angles  of  10°,  30°,  50°,  70°,  and  90°. 

Unnormalized  values  of  Sn(180°)  are  shown  in  Figure  5.12;  in  these  and 
subsequent  results  the  local  maximum  in  the  backscattering  signal  for  10°  is 
about  twice  that  for  the  other  angles.  In  succeeding  figures  the  local  maximum 
has  been  normalized  to  1.0  for  comparing  relative  changes  in  the 

O 

backscattering  signal.  Normalized  values  of  Su(180  )  are  shown  in  figure  5.13. 
For  the  zenith  angles  30°  and  90°  the  relative  difference  between  the  local 
maximum  and  minimum  is  small,  and  it  is  largest  for  50°  and  70°.  The  shapes 
of  the  curves  are  somewhat  different  for  values  of  Shh(180°)  and  Sw(180°) 
(shown  in  Figures  5.14  and  5.15).  The  relative  difference  between  the  local 
maximum  and  minimum  for  the  5^(180°)  curves  are  more  pronounced  for  the 
smallest  zenith  angles.  For  Sw(180°)  the  relative  difference  between  local 
maximum  and  minimum  is  similar  for  all  angle  except  90°,  but  the  noticeable 
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Figure  5.12  Backscattering  of  unpolarized  light  by  ice  columns  of  increasing 
size  parameter  at  various  zenith  angles.  The  short-dashed  line 

( - )  represents  a  zenith  angle  of  10°,  the  dot-dash  line  ( . ) 

zenith  angle  of  30  ,  the  dotted  line  a  zenith  angle  of  50  ,  the 
solid  line  a  zenith  angle  of  70  ,  and  the  long-dashed  line  ( — 
- )  a  zenith  angle  of  90°. 
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Figure  5.13  Backscattering  of  unpolarized  light  by  ice  columns  of  increasing 
size  parameter  at  various  zenith  angles.  Same  as  figure  5.10 
except  the  local  maximum  value  in  the  curves  have  been 

normalized  to  1.0.  The  short-dashed  line  ( - )  represents  a 

zenith  angle  of  10°,  the  dot-dash  line  ( . )  a  zenith  angle  of 

30  ,  the  dotted  line  a  zenith  angle  of  50°,  the  solid  line  a  zenith 

angle  of  70  ,  and  the  long-dashed  line  ( - )  a  zenith  angle 

of  90°. 
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Figure  5.14  Backscattering  of  horizontally  polarized  light  by  ice  columns  of 
increasing  size  parameter  at  various  zenith  angles.  The  short- 
dashed  line  (-  -  -)  represents  a  zenith  angle  of  10  ,  the  dot-dash 

line  ( . )  a  zenith  angle  of  30  ,  the  dotted  line  a  zenith  angle 

of  50°,  the  solid  line  a  zenith  angle  of  70  ,  and  the  long-dashed 
line  ( - )  a  zenith  angle  of  90  . 
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Figure  5.15  Backscattering  of  vertically  polarized  light  by  ice  columns  of 
increasing  size  parameter  at  various  zenith  angles.  The  short- 
dashed  line  (-  -  -)  represents  a  zenith  angle  of  10  ,  the  dot-dash 

line  ( . )  a  zenith  angle  of  30°,  the  dotted  line  a  zenith  angle 

of  50  ,  the  solid  line  a  zenith  angle  of  70°,  and  the  long-dashed 
line  ( - )  a  zenith  angle  of  90  . 
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Figure  5.16  Backscattering  by  ice  columns  of  increasing  size  parameter  at 
zenith  angle  of  70  .  The  dashed-line  represents  unpolarized 
light,  the  dotted  line  represents  horizontally  polarized  light,  and 
the  solid  line  represents  vertically  polarized  light. 
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feature  is  the  narrowing  of  the  local  maximum  feature  at  the  two  largest  zenith 
angles.  A  comparison  of  the  backscattering  signal  at  70°  for  all  three  states  of 
polarization  is  shown  in  Figure  5.16.  The  most  pronounced  local  minimum  is 
calculated  to  appear  near  a  size  parameter  of  2.0  for  Sw(180°).  The 
comparison  of  the  backscattering  signals  in  Figures  5.12  through  5.16  agrees 
with  one’s  intuition  that  linearly  polarized  microwave  radar  (in  particular  Sw) 
would  produce  better  definition  in  the  backscattering  curves.  This  is  also 
evident  in  backscattering  by  ice  plates. 

Figures  5.17  through  5.20  show  the  dependence  of  backscattering  by  the 
hexagonal  plate  (aspect  ratio  of  3.2:1)  on  size  parameter.  The  size  parameter 
calculation  is  again  based  on  an  equal-volume  sphere.  The  values  for  size 
parameter  range  from  0.6  to  2.4  and  correspond  to  plate  diameters  of  0.8  to 
3.3  mm.  The  five  lines  in  the  figures  represent  zenith  angles  of  10°,  30°,  50°, 
70°,  and  90°.  Normalized  values  of  Sn(180°)  are  shown  in  Figure  5.17. 

Unlike  the  results  for  the  hexagonal  column,  the  backscattering  minimum  for 
the  plate  is  spread  over  a  larger  range  of  size  parameters,  and  it  occurs  at  a 
lower  size  parameter  for  higher  zenith  angles.  The  relative  difference  between 
the  local  maximum  and  local  minimum  also  decreases  with  increasing  zenith 
angle.  For  values  of  Shh(180°)  and  Sw(180°)  (shown  in  Figure  5.18  and  5.19) 
the  shapes  of  the  curves  are  again  dissimilar.  The  main  difference  between 
Shh(180o)  and  Su(180°)  is  a  slight  deepening  of  the  backscattering  minimum  at 
a  zenith  angle  of  50°.  The  noticeable  differences  between  Sw(180°)  and 
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Figure  5.17  Backscattering  of  unpolarized  light  by  ice  plates  of  increasing  size 
parameter  at  various  zenith  angles.  The  short-dashed  line  (-  -  -) 
represents  a  zenith  angle  of  10  ,  the  dot-dash  line  (-•-•-)  a 
zenith  angle  of  30  ,  the  dotted  line  a  zenith  angle  of  50  ,  the 
solid  line  a  zenith  angle  of  70  ,  and  the  long-dashed  line  ( — 
- )  a  zenith  angle  of  90°. 
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Figure  5.18  Backscattering  of  horizontally  polarized  light  by  ice  plates  of 
increasing  size  parameter  at  various  zenith  angles.  The  short- 
dashed  line  (-  -  -)  represents  a  zenith  angle  of  10  ,  the  dot-dash 

line  ( . )  a  zenith  angle  of  30°,  the  dotted  line  a  zenith  angle 

of  50°,  the  solid  line  a  zenith  angle  of  70°,  and  the  long-dashed 
line  ( - )  a  zenith  angle  of  90°. 
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Figure  5.19  Backs  cattering  of  vertically  polarized  light  by  ice  plates  of 

increasing  size  parameter  at  various  zenith  angles.  The  short- 
dashed  line  (-  -  -)  represents  a  zenith  angle  of  10°,  the  dot-dash 

line  ( . )  a  zenith  angle  of  30°,  the  dotted  line  a  zenith  angle 

of  50°,  the  solid  line  a  zenith  angle  of  70°,  and  the  long-dashed 
line  ( - )  a  zenith  angle  of  90°. 
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Figure  520  Backscattering  by  ice  plates  of  increasing  size  parameter  at  zenith 
angle  of  50°.  The  dashed-line  represents  unpolarized  light,  the 
dotted  line  represents  horizontally  polarized  light,  and  the  solid 
line  represents  vertically  polarized  light. 
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Su(180°)  are  the  shallower  minimum  at  zenith  angle  of  50°  and  the  presence 
of  the  second  backscattering  minimum  for  all  zenith  angles  except  10°  and  90° 

(a  shallow  minimum  occurs  at  70°,  but  it  is  out  of  range  of  the  figure).  A 
comparison  of  the  backscattering  signal  at  50°  for  all  three  states  of 
polarization  is  shown  in  Figure  5.20.  In  contrast  with  the  results  for  the 
hexagonal  column,  SHH(180°)  produces  the  largest  relative  difference  between 
the  local  maximum  and  minimum. 

The  backscattering  minima  for  ice  crystals  have  an  added  feature  over 
their  counterparts  for  spherical  raindrops:  a  dependence  on  zenith  angle.  This 
dependence  disappears  if  the  ice  crystals  tumble  as  they  fall;  nonetheless,  this 
could  provide  information  about  turbulence  in  the  vicinity  of  the  crystals. 

Thus,  backscattered  signals  collected  while  the  radar  sweeps  through  zenith 

angles  may  be  useful  In  remotely  sensing  the  atmospheric  environment. 

timra  1  firm  wiiiifn  Out  mlm  of  /nmm)1o§y  would  fmeflt  by  Iwvbig 
a  better  understanding  of  the  crystal  characteristics  within  ice  clouds:  climate 
modeling,  atmospheric  chemistry,  and  cloud  physics  (Vogelmann,  et  al.,  1990). 
The  transfer  of  infrared  radiation  through  cirrus  clouds  is  an  issue  in  the 
debate  over  global  warming.  Modeling  scattering  by  cirrus  clouds  can  be 
improved  if  the  size  and  shape  distribution  of  the  ice  crystals  are  known.  The 
scavenging  or  removal  of  atmospheric  aerosols,  including  sulfates  that 
contribute  to  acid  precipitation,  depends  on  the  shape  of  the  crystal.  The 
scavenging  efficiency  of  ice  crystals  appears  to  increase  with  surface  to  volume 
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ratio  (Dungey,  1976).  Crystal  growth  and  aggregation  are  important 
mechanisms  in  the  precipitation  process.  By  analyzing  backscattering  signals 
from  different  levels  within  the  atmosphere  the  growth  of  the  ice  crystals  may 
be  estimated. 

Although  identification  of  size  and  shape  distributions  of  ice  crystals 
using  this  procedure  seems  feasible,  several  contributing  technological  factors 
are  still  uncertain.  The  technology  required  to  build  millimeter-wavelength 
radar  is  relatively  new.  Present  research  and  development  of  millimeter  wave 
tubes  is  expected  to  result  in  construction  of  more  powerful  transmitters 
(Lhermitte,  1990),  thereby  increasing  the  sensitivity  of  the  radar.  The  size 
range  of  ice  crystals  in  this  chapter  is  typical  of  falling  snow.  To  identify  the 
size  distribution  of  smaller  ice  crystals,  a  shorter  wavelength  is  required.  A 
recently  developed  1.4  mm  (215  GHz)  wavelength  radar  (Mead  et  al.,  1989) 
has  the  potential  for  observing  the  first  minimum  in  backscattering  by  smaller 
ice  crystals;  however,  at  this  wavelength  attenuation  by  atmospheric  water 
vapor  may  (depending  of  the  transmission  power)  limit  the  system’s 
performance. 

Identification  of  ice  crystal  size  distribution  depends  on  several  other 
factors.  First,  the  size  distribution  is  assumed  to  be  continuous.  A 
discontinuous  size  distribution  would  produce  gaps  in  the  Doppler  shift 
spectrum  which  could  be  mistaken  as  a  backscattering  minimum.  Second, 
crystal  aggregation  should  be  a  minimum.  Backscattering  calculations  are 
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made  only  for  plates  and  columns;  scattering  by  aggregates  may  irretrievably 
obscure  the  desired  backscattering  minima.  Third,  the  optical  depth  of  the  ice 
crystals  is  not  too  large,  otherwise  the  backscattered  signal  will  be  a  result  of 
multiple  scattering  events.  The  procedure  for  identifying  ice  crystals  is  based 
on  single  scattering. 

In  this  chapter,  the  ability  of  the  coupled-dipole  to  calculate 
backscattering  by  nonspherical  particles  was  demonstrated.  For  spheres 
smaller  than  the  size  which  produces  the  first  backscattering  minimum,  small 
shape  variations  in  the  dipolar  array  do  not  affect  the  accuracy  of  the 
backscattering  calculation.  Conversely,  when  the  spheres  are  larger  than  the 
size  which  produces  the  first  backscattering  minimum,  backscattering 
calculations  lose  accuracy.  For  particles  with  a  large  aspect  ratio  (at  least  5:1) 
and  small  size  parameter,  the  backscattering  calculation  is  relatively  insensitive 
to  the  exact  shape  of  the  particle,  i.e.  hexagonal  column,  cylinder,  etc.  Finally, 
backscattering  calculations  for  hexagonal  columns  and  plates  have  been 
presented.  Although  this  is  not  a  definitive  work  in  radar  meteorology,  it  is 
feasible  that  a  size  distribution  of  ice  crystals  could  be  estimated  using  remote 
sensing  techniques.  This  can  also  be  scaled  down  for  observing  cirrus  crystals 
with  an  instrument  operating  at  a  shorter  wavelength.  In  general,  the  size  of 
the  ice  crystals  will  dictate  the  wavelength  that  will  be  best  suited  for  obtaining 
the  first  backscattering  minimum. 
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Chapter  6 


SUMMARY  AND  CONCLUSION 


The  problem  of  light  scattering  by  arbitrary  particles  can  be  approached 
a  number  of  ways.  One  option  too  often  taken  is  to  assume  that  the  particles 
are  spherical.  This  assumption  simplifies  the  computational  task  by  allowing 
Mie  theory  to  be  used  for  the  calculations;  however,  this  oversimplification 
may  lead  to  unacceptable  errors.  The  other  option  is  to  appreciate  the  shape 
of  the  particles  and  choose  an  appropriate  computer  model  to  calculate  the 
particle’s  scattering  properties. 

As  described  in  this  dissertation,  the  coupled-dipole  method  is 
applicable  for  modeling  a  considerable  variety  of  arbitrary  particles.  First,  a 
particle  of  any  shape  can  be  modeled  using  this  method.  Second,  several 
solution  techniques  are  available  depending  on  the  physical  characteristics  of 
the  particle.  If  the  size  parameter  of  the  particle  is  small  ( <  2  or  3  depending 
on  shape  or  refractive  index)  or  if  the  particle  is  to  be  modeled  in  numerous 
orientations,  the  matrix  inversion  technique  may  be  used.  If  the  relative 
refractive  index  of  the  particle  is  not  too  large  ( <  2),  computer  memory 
limitations  can  be  avoided  by  using  the  scattering-order  technique. 

Another  consideration  in  the  coupled-dipole  method  is  the  relationship 
between  the  polarizability  of  the  discrete  dipolar  subunits  and  the  relative 


95 


refractive  index  of  the  particle  being  modeled.  In  most  papers  the  Clausius- 
Mosotti  (CM)  relation  or  the  CM  relation  with  a  radiative  reaction  term  has 
been  used.  Since  the  formulation  of  the  coupled-dipole  method  is  based  on 
subunits  that  radiate  as  dipoles,  the  conversion  of  relative  refractive  index  to 
electric  dipole  polarizability  (using  the  first  term  from  Mie  theory  referred  to 
as  Doyle’s  method)  was  incorporated  into  the  method.  The  results  using 
Doyle’s  method  show  a  better  agreement  with  Mie  theory  than  using  the  other 
schemes.  Further  calculations  showed  that  Doyle’s  method  can  be  used  near  a 
resonance  mode  (Frohlich  mode)  for  small  particles  and  that  the  results 
compared  favorable  with  measured  data. 

Input  parameters  for  the  coupled-dipole  method  have  been  added  to 
accommodate  a  variety  particle  shapes  and  sizes.  Subroutines  are  available  for 
constructing  arrays  of  dipolar  units  that  represent  spheres,  spheroids, 
rectangular  solids,  cylinders,  and  hexagonal  crystals.  The  capability  of 
specifying  the  size  of  particle  has  also  been  included.  Output  provides  a 
number  of  useful  scattering  parameters:  Mueller  matrix  elements  for  any 
scattering  angle  as  well  as  scattering,  absorption,  and  extinction  cross  sections. 

The  ability  of  the  coupled-dipole  method  for  determining  backscattering 
calculations  was  investigated.  A  simple  model  was  presented  to  show  why 
scattering  in  the  backward  direction  is  most  sensitive  to  particle  size  and  shape. 
The  sensitivity  was  shown  to  be  a  result  of  the  phase  relation  of  scattered 
waves  from  the  dipolar  subunits.  Backscattering  values  using  the  coupled- 
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dipole  method  compare  favorably  with  Mie  theory  for  spheres  which  are 
smaller  than  the  size  parameter  associated  with  the  first  backscattering 
minimum  (e.g.,  size  parameter  of  1.6  for  relative  refractive  index  of 
1.33  +  i0.05).  At  this  and  larger  size  parameters,  backscattering  becomes 
extremely  shape  dependent,  and  modeling  results  from  the  coupled-dipole 
method  lose  precision  because  an  array  of  dipolar  subunits  arranged  on  a 
simple  cubic  lattice  cannot  adequately  represent  a  solid  sphere.  Small 
variations  in  shape  were  shown  to  strongly  affect  backscattering  calculations 
for  spheres  with  size  parameters  which  are  near  the  size  parameter  associated 
with  the  first  backscattering  minimum.  Agreement  between  the  coupled-dipole 
method  and  Mie  theory  can  be  improved  by  increasing  the  number  of  dipolar 
subunits  in  the  array;  however,  this  also  increases  computer  memory 
requirements  and  computation  time. 

Backscattering  by  equal-volume  particles  was  investigated  using  the 
coupled-dipole  method.  A  cylinder,  rectangular  solid,  prolate  spheroid,  and 
hexagonal  column  of  equal  volume  and  with  aspect  ratio  5:1  were  modeled, 
and  the  results  were  compared  with  those  of  an  equal-volume  sphere.  Despite 
the  difference  in  cross-sectional  shape,  the  elongated  particles  displayed  similar 
scattering  properties  for  three  different  refractive  indices.  However,  the  size 
parameter  of  the  particles  was  only  1.0.  If  the  particles  were  larger  than  the 
size  which  corresponds  to  the  first  backscattering  minimum,  backscattering  by 
the  elongated  particles  might  not  agree  so  well  because  of  the  extreme  shape 
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dependence  as  discussed  above.  The  average  backscattering  value  for  the 
elongated  shapes  (this  simulated  random  orientation  of  the  particles)  was 
approximately  equal  to  that  for  the  equal-volume  sphere;  however,  agreement 
decreased  with  increasing  refractive  index. 

Finally,  backscattering  of  94  GHz  radar  by  ice  crystals  was  examined. 
The  accuracy  of  the  coupled-dipole  method  in  this  application  was  shown  to  be 
good  by  comparing  calculations  of  backscattering  by  ice  spheres  with  Mie 
theory.  Backscattering  by  randomly  oriented  ice  columns  and  plates  was 
compared  with  backscattering  by  equal-volume  ice  spheres.  For  small  size 
parameters  the  backscattering  values  are  in  fair  agreement;  however,  for  larger 
size  parameters  they  differ  by  as  much  as  a  factor  of  six. 

The  orientation  of  the  ice  columns  and  plates  was  then  made  more 
restrictive  by  modeling  them  as  if  they  fell  through  the  atmosphere  without 
tumbling.  The  resulting  backscattering  calculations  now  contained  a 
dependence  on  zenith  angle.  The  first  backscattering  minima  are  present  for 
both  type  crystals,  but  with  noticeable  differences.  The  relative  difference 
between  the  local  maximum  and  the  first  backscattering  minimum  for  the  ice 
plate  is  larger  than  that  for  the  column,  and  the  backscattering  minimum  for 
the  plate  is  spread  out  over  a  larger  range  of  size  parameters.  The  use  of 
horizontally  and  vertically  polarized  waves  also  enhance  the  magnitude  of  the 
first  backscattering  minima. 
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Backscattering  minima  by  raindrops  can  be  detected,  and  the 
information  is  used  to  determine  additional  features  of  the  atmosphere.  These 
variations  in  backscattering  minima  by  ice  crystals  at  94  GHz  may  be  used  to 
remotely  estimate  their  size  and  shape  distributions.  Knowledge  of  the  ice 
crystal  characteristics  would  benefit  the  areas  of  climate  modeling,  atmospheric 
chemistry,  and  cloud  physics.  Present  research  and  development  of  millimeter 
radar  technology  will  provide  increased  sensitivity,  which  makes  this  scheme 
for  identifying  ice  crystals  even  more  promising  (Lhermitte,  1988). 

In  conclusion,  scattering  by  arbitrary  particles  can  be  approximated  by 
assuming  their  scattering  properties  are  similar  to  that  of  equal-volume 
spheres.  This  practice  has  advantages  and  disadvantages.  The  calculations  for 
scattering  by  equal-volume  spheres  are  quick  and  inexpensive,  and  in  many 
instances  they  provide  a  reasonable  first-order  estimate  of  the  scattering 
properties  of  the  arbitrary  particles  that  may  not  have  been  initially  apparent. 
But  by  using  the  equal-volume  sphere  approach,  additional  scattering 
information  may  be  overlooked. 

It  was  shown  that  scattering  by  an  equal-volume  sphere  agrees  favorably 
with  that  for  small,  tumbling  ice  crystals  at  94  GHz  (size  parameter  less  than 
0.8).  Thus  in  this  case  the  equal-volume  sphere  approximation  will  save 
computer  time  and  costs  while  providing  reasonable  results.  However,  for  size 
parameters  larger  than  1.0,  the  equal-volume  sphere  approximation  will  lead  to 
large  errors.  Moreover,  when  the  ice  crystals  are  assumed  to  be  not  tumbling, 
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the  scattering  results  take  on  additional  features.  The  backscattered  signal 
becomes  dependent  on  zenith  angle  and  on  the  polarization  of  the  incident 
wave.  These  properties  of  scattering  may  help  to  identify  the  shape  and  size 
distributions  of  ice  crystals,  but  they  will  never  be  present  for  spheres. 
Therefore,  while  calculating  scattering  by  spheres  may  be  a  simple  task  and 
sometimes  a  fairly  accurate  approximation,  important  information  could  be 
omitted. 

The  next  problem  to  be  solved  is  that  of  inverting  the  backscattered 
signal  to  useful  information  about  the  size  and  shape  of  the  ice  crystals. 
Variations  in  backscattering  for  ice  crystals  of  different  aspect  ratios  should 
also  be  examined.  The  refractive  index  of  ice  at  94  GHz  is  sufficiently  large 
to  induce  higher  order  multipoles  in  the  ice  crystals.  It  may  be  advantageous 
to  incorporate  the  calculation  of  these  multipoles  in  the  coupled-dipole  method 
to  improve  the  accuracy  of  scattering  calculations. 
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Appendix  A 


This  appendix  has  been  accepted  for  publication  in  the  Journal  of  the 
Optical  Society  of  America  A.  A  slightly  modified  version  of  the  paper 
is  given  here. 

LIGHT  SCATTERING  BY  NONSPHERICAL  PARTICLES: 

A  REFINEMENT  TO  THE  COUPLED-DIPOLE  METHOD 


Clifton  E.  Dungeyt  and  Craig  F.  Bohren 
Department  of  Meteorology,  The  Pennsylvania  State  University 
University  Park,  Pennsylvania  16802 


Abstract 

In  the  coupled-dipole  method  an  arbitrary  particle  is  modeled  as  an 
array  of  N  polarizable  subunits  each  of  which  gives  rise  to  only  electric  dipole 
radiation.  The  Clausius-Mosotti  relation  is  widely  used  to  calculate  the 
polarizability  of  the  subunits  that  corresponds  to  the  dielectric  function  of  the 
particle  that  the  array  represents.  In  this  paper  we  replace  the  Clausius- 
Mosotti  relation  with  an  exact  expression  for  the  electric  dipole  polarizability 
and  find  improvement  in  extinction  calculations  for  spheres  as  compared  with 
Mie  theory.  Near  a  Frohlich  frequency  the  coupled-dipole  method  yields 
extinction  cross  sections  for  spheres  and  spheroids  that  compare  favorably  with 
the  method  of  continuous  distribution  of  ellipsoids  and  measured  values. 


+Captain  Dungey  is  assigned  to  Penn  State  through  the  graduate  meteorology 
program  of  the  Air  Force  Institute  of  Technology. 
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Introduction 

The  coupled-dipole  method  was  apparently  first  applied  by  Purcell  and 
Pennypacker  (1973)  to  calculate  approximate  cross  sections  of  nonspherical 
particles.  Partly  owing  to  the  enormous  computer  storage  and  CPU  time 
required  for  modeling  particles  with  large  size  parameters,  this  method  never 
has  had  a  great  following.  But  in  the  past  few  years,  with  faster  computers 
and  more  efficient  programming  techniques,  an  increasing  number  of  people 
have  begun  using  this  relatively  simple  technique.  For  example,  Draine  (1988) 
has  used  it  to  study  scattering  by  interstellar  grains;  Goedecke  and  O’Brien 
(1988)  and  Flatau  et  al.  (1990)  have  examined  scattering  by  ice  crystals; 

S.  B.  Singham  et  al.  (1986a)  calculated  differential  scattering  by  chiral 
particles;  and  Varadan  et  al.  (1989)  computed  scattering  by  anisotropic 
particles. 

The  coupled-dipole  method  has  remained  essentially  unchanged  since  its 
inception.  An  arbitrary  particle  is  divided  into  an  array  of  N  subunits  located 
on  a  simple  cubic  lattice  (Figure  A.l).  The  dipolar  subunits  are  sufficiently 
small  to  give  rise  to  only  electric  dipole  radiation.  Total  scattering  is  then 
calculated  by  summing  the  waves  scattered  by  each  dipolar  subunit  excited  by 
the  incident  wave  and  by  the  waves  scattered  to  it  from  all  its  neighbors.  The 
coupled-dipole  method  originally  was  formulated  (Purcell  and  Pennypacker, 
1973)  heuristically;  however,  a  more  formal  mathematical  derivation  and  a 
short  review  of  the  method  has  been  published  recently  (Lakhtakia,  1990). 

Some  recent  modifications  made  to  the  coupled-dipole  method  include 
more  efficient  means  of  calculating  scattering  by  randomly  oriented  particles 
(Singham  et  al.,  1986b)  and  improved  solution  algorithms  (Draine,  1988;  Flatau 
et  al.,  1990;  Chiapetta,  1980;  Singham  and  Bohren,  1988)  that  permit  larger 
values  of  N.  For  scattering  calculations  in  this  paper  we  rely  on  the  matrix 
inversion  technique  (Singham  and  Salzman,  1986)  to  solve  the  interaction 


Figure  A.1  A  spherical  particle  which  has  been  represented  by  an  array  of 
136  dipolar  subunits.  The  effective  radius  of  the  sphere  is 
determined  by  a*  =  (3N/4;r)w,  where  the  dipolar  subunits  are 
assumed  to  have  unit  volume. 
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matrix  directly.  Although  N  is  limited  to  300  in  this  technique  because  of  the 
need  to  store  the  matrix,  we  are  guaranteed  a  solution  provided  that  the  matrix 
is  algorithmically  nonsingular.  The  matrix  inversion  method  also  permits  the 
efficient  calculation  of  extinction  cross  sections  of  nonspherical  particles  in 
random  orientations  (Singham  et  al.,  1986b).  In  this  paper  we  treat  only 
isotropic  scatterers. 

An  integral  part  of  the  coupled-dipole  method  is  an  effective-medium 
theory  that  provides  the  relationship  between  the  polarizability  of  the 
individual  dipolar  subunits  and  the  bulk  dielectric  function  of  the  continuous 
medium  that  the  array  represents.  The  most  widely  used  effective-medium 
theory  for  determining  the  polarizability  has  been  the  Clausius-Mosotti  (CM) 
relation,  in  some  instances  with  slight  modifications  (e.g.  Draine,  1988). 

In  next  section,  following  Doyle  (1989),  we  calculate  the  electric  dipole 
polarizability  using  an  exact  expression.  Using  this  expression  in  the  coupled- 
dipole  method  generally  leads  to  calculated  extinction  cross  sections  more  in 
accord  with  Mie  theory. 

In  final  section  we  test  Doyle’s  method  on  particles  in  the  neighborhood 
of  a  Frohlich  frequency.  This  is  a  strong  absorption  mode  highly  dependent  on 
both  particle  size  and  shape.  Extinction  cross  sections  near  the  Frohlich  mode 
are  computed  for  spheres  and  spheroids  by  using  the  coupled-dipole  method 
and  are  compared  with  results  from  Mie  theory,  the  method  of  continuous 
distribution  of  ellipsoids  (CDE),  and  measured  values  (Huffman  and  Bohren, 
1980).  The  nature  of  the  shape-dependent  extinction  determined  by  the 
coupled-dipole  method  is  used  to  analyze  small  discrepancies  between  results 
from  modeling  spheres  with  the  coupled-dipole  method  and  Mie  theory. 
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Effective-Medium  Theory 

When  dealing  with  a  composite  material  consisting  of  small  grains 
embedded  in  a  homogeneous  matrix,  one  must  measure  or  estimate  the 
dielectric  function  from  a  theory.  If  the  optical  properties  of  the  embedded 
grains  and  the  matrix  are  known,  one  uses  an  effective-medium  theory  to 
estimate  the  bulk  characteristics  of  the  composite.  This  problem  is  the 
converse  of  that  faced  by  those  using  the  coupled-dipole  method  for  which  the 
bulk  dielectric  functions  of  the  particle  (composite)  and  surrounding  medium 
are  known  and  the  optical  properties  of  the  dipolar  subunits  (grains)  must  be 
calculated. 

The  Maxwell  Garnett  relation  is  an  effective-medium  theory  that  can  be 
used  to  determine  the  dielectric  function  of  the  dipolar  subunits.  By  assuming 
spherical  dipolar  subunits,  we  can  reduce  the  Maxwell-Gamett  relation 
(Bohren  and  Huffman,  1983,  Sect  8.5): 

2(1  -/)e2„  -  (2+/)e„e  (A.1) 

(l-/)e  -  (l+2/)enl 

where  edip,  e,  and  em  are  the  dielectric  functions  of  the  dipolar  subunits,  the 
bulk  particle,  and  its  surrounding  medium;  and  f  is  the  volume  filling  factor. 
For  spherical  dipolar  subunits  on  a  simple  cubic  lattice  with  diameters  equal  to 
the  lattice  spacing  f  =  n/6. 

Purcell  and  Pennypacker  (1973)  used  the  CM  relation  to  obtain  the 
complex  polarizability  a  of  the  dipolar  subunits  from  the  dielectric  function  of 
the  material: 
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a  _  3(g  ~  gm)  (A.2) 

N(e  +  2em) 

where  N  is  the  number  of  dipolar  subunits  per  unit  volume.  Although  the 
Maxwell-Gamett  and  CM  relations  yield  different  physical  parameters,  their 
formal  derivations  have  been  shown  to  be  essentially  identical  (Barker,  1973). 
The  relationship  between  a  and  edip  is  obtained  by  combining  equations  (1) 
and  (2)  producing  another  form  of  the  CM  relation: 


a  -  47t£3 


(*Jip 

(Bdip 


g  J 
20 


(A.3) 


where  a  is  the  radius  of  the  dipolar  subunits.  Equations  (1)  and  (3)  will  be 
referred  to  later,  but  next  we  briefly  examine  a  shortcoming  in  the  CM  relation 
when  it  is  used  in  the  coupled-dipole  method. 

The  optical  theorem  states  that  the  extinction  cross  section  Cext  for  small 
spheres  is  determined  from  the  imaginary  part  of  its  polarizability.  When  e  is 
purely  real  the  CM  relation  produces  a  purely  real  a.  Since  Cext  cannot  be 
zero  because  the  incident  wave  must  be  attenuated  by  scattering,  a  must  be 
complex.  Draine  (1988)  and  Goedecke  and  O’Brien  (1988)  provide  methods 
that  satisfy  this  criterion,  we  introduce  a  third. 

When  investigating  the  optical  properties  of  small  metal  spheres 
suspended  in  a  transparent  medium,  Doyle  (1989)  used  an  exact  expression  for 
the  electric  dipole  polarizability 


a , 


(A.4) 


where  k  =  2jta./\  and  a,  is  the  electric  dipole  coefficient  from  Mie  theory: 
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a 


(A.5) 


with  v>i  and  £x  being  Riccati-Bessel  functions  and  x  =  ka.  In  Doyle’s 
application  m  is  the  ratio  of  the  complex  index  of  refraction  of  the  metal 
spheres  to  that  of  the  medium;  for  the  coupled-dipole  method  it  is  the  ratio  of 
the  index  of  refraction  of  the  dipolar  subunits  to  that  of  the  medium.  When 
using  this  exact  expression  for  the  electric  dipole  term,  Doyle  calculated 
reflectances  of  a  composite  that  agreed  better  with  experimental  values  than  if 
only  the  first  term  in  the  series  had  been  used. 

The  coupled-dipole  method  operates  on  the  principle  that  the  dipolar 
subunits  are  sufficiently  small  that  they  give  rise  to  only  electric  dipole 
radiation.  It  was  therefore  appropriate  to  incorporate  equations  A.4  and  A.5 
into  the  coupled-dipole  method.  In  what  we  now  refer  to  as  the  Doyle 
expression,  the  polarizability  of  the  dipolar  subunits  are  calculated  using  this 
exact  expression  for  the  electric  dipole  polarizability.  The  value  for  m  in 
equation  A.5  is  obtained  from  equation  A.l  where  m2  =  edip/em.  In  the  small- 
particle  limit,  the  Doyle  expression  reduces  to  the  CM  relation. 

To  guarantee  a  complex  a,  Draine  used  a  radiative  reaction  term  that  is 
contained  in  the  Doyle  expression.  If  ax  were  expanded  in  a  series  of  the  size 
parameter  x,  the  radiative  reaction  term  would  be  the  third  in  the  series.  The 
Doyle  expression  and  Draine’s  radiative-reaction  factor  are  similar  in  that  they 
modify  the  entire  interaction  matrix.  Goedecke  and  O’Brien  (1988)  use  a  self¬ 
term  correction  that  modifies  only  the  diagonal  elements  of  the  matrix.  We 
now  compare  results  from  the  coupled-dipole  method  by  using  the  CM 
relation,  Draine’s  radiation  correction,  and  Doyle’s  expression.  In  our  model 
all  the  dipolar  subunits  are  spherical;  most  of  the  particles  (dipolar  arrays)  that 
are  modeled  are  also  spherical.  In  our  discussions  we  try  to  distinguish 
between  the  two. 
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In  this  section  dipolar  arrays  representing  spheres  are  modeled  and  the 
results  are  compared  with  BHMIE  (Bohren  and  Huffman,  1983,  Appendix  A). 
For  the  first  comparison  the  arrays  range  in  size  from  27  to  251  dipolar 
subunits.  If  the  spacing  of  the  lattice  is  considered  to  be  1  du  (dipole  unit) 
then  the  corresponding  effective  radii  for  the  modeled  spheres  range  from  1.86 
to  3.91  du.  A  wavelength  of  12.5  du  was  used.  Figure  A.2  shows  relative 
values  of  for  two  refractive  indices  using  the  three  schemes  to  determine 
the  dipolar  subunit  polarizability.  For  this  comparison  the  Doyle  expression 
agrees  best  with  Mie  theory. 

One  input  variable  remained  constant  for  this  comparison:  the  size 
parameter  of  the  dipolar  subunits.  Purcell  and  Pennypacker  (1973)  state  that  if 
ka  <  0.35,  the  coarseness  of  the  array  would  not  cause  serious  errors  in  the 
calculations.  Yung  (1978)  claims  reliable  results  when  ka  <  0.17.  By  also 
considering  the  wavelength  of  the  electric  field  within  the  dipolar  subunit 
Draine’s  guidance  (1988)  suggests  ka|m[  <  0.5  when  an  accuracy  of  10%  is 
desired,  where  |m|  is  the  modulus  of  the  complex  refractive  index;  for  this 
comparison  these  values  were  0.33  and  0.43. 

In  the  second  comparison  the  size  of  the  dipolar  subunit  was  varied. 

The  range  of  ka|m|  was  from  0.13  to  1.07.  The  sphere  was  represented  by  an 
array  of  136  dipolar  subunits;  the  wavelength  was  12.5  du.  As  seen  from 
Figure  A.3,  when  the  size  parameter  of  the  subunits  is  small,  all  three  schemes 
for  determining  polarizability  compare  well  with  Mie  theory.  Over  the  entire 
range  the  Doyle  expression  is  the  most  accurate.  Note  the  increasing 
overprediction  by  all  three  schemes  above  ka|m|  =  0.8,  which  occurs  as  Qert  of 
the  bulk  sphere  is  reaching  a  maximum  value  of  3.7.  Not  only  does  the  Doyle 
expression  overpredict  least,  it  is  the  only  one  to  respond  to  the  falling  values 
of  Qext  at  the  largest  size  parameter  indicated. 
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Size  Parameter  of  Sphere 


Figure  A.2  Comparison  between  Mie  theory  and  the  coupled-dipole  method 
using  different  schemes  to  calculate  polarizability  while  varying 
the  number  of  dipolar  subunits.  Calculations  using  the  Doyle 

expression  are  represented  by  ( - ),  Draine’s  radiative 

reaction  term  by  ( •  •  •  • ),  and  the  Clausius-Mosotti  relation  by 

( - ).  n  varies  from  27  to  251;  size  of  dipolar  subunit  remains 

constant.  Results  for  two  refractive  indices  are  shown. 
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Size  Parameter  of  Dipolar  Subunit 


Figure  A3  Comparison  between  Mie  theory  and  the  coupled-dipole  method 
using  different  schemes  to  calculate  polarizability  while  varying 
the  size  of  the  dipolar  subunits.  Calculations  using  the  Doyle 

expression  are  represented  by  ( - ),  Draine’s  radiative 

reaction  term  by  ( •  *  •  • ),  and  the  Clausius-Mosotti  relation  by 
( - ).  N  remains  constant  (136  dipolar  subunits). 
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The  index  of  refraction  was  varied  in  the  third  comparison.  An  array  of 
136  dipolar  subunits  was  modeled;  the  effective  radius  of  the  particle  remained 
constant  and  the  wavelength  was  12.5  du.  Nine  refractive  indices  were 
selected  from  earlier  papers  (Purcell  and  Pennypacker,  1973;  Draine,  1988) 
and  from  various  practical  applications.  As  seen  in  Table  A.1,  agreement 
between  the  coupled-dipole  method  and  Mie  theory  generally  varies  inversely 
with  ka  |m| .  Among  the  different  schemes  the  Doyle  expression  gives  the  best 
overall  results.  For  instance,  no  improvement  is  noted  for  m  =  1.44  +  i0.26 
but  the  coupled-dipole  method  already  compares  to  within  1%  of  Mie  theory. 
At  the  two  largest  refractive  indices  the  Doyle  expression  does  not  compare 
best  with  Mie  theory.  It  may  be  that  the  other  schemes  are  near  a  crossover 
point  as  was  the  case  in  Figure  A.2  near  ka|mj  =  0.8.  However,  the  use  of  the 
coupled-dipole  method  in  this  particular  application  is  not  recommended 
because  the  higher  order  multipoles  associated  with  the  large  refractive  index 
are  undoubtedly  introducing  error  in  the  calculations.  A  better  comparison 
would  be  made  with  a  much  array  with  a  larger  number  of  dipoles; 
unfortunately,  the  array  size  is  limited  in  the  matrix-inversion  method. 

When  the  coupled-dipole  method  is  used  appropriately  (ka|m|  <  0.5) 
the  Doyle  expression  yields  extinction  efficiencies  that  agree  better  with  Mie 
theory  than  the  other  two  schemes  do.  Incorporation  of  the  Doyle  expression 
into  the  coupled-dipole  method  adds  negligible  computation  time. 

The  refractive  indices  used  in  the  previous  comparisons  are  appropriate 
only  for  regions  in  which  surface  modes  are  not  excited.  A  more  stringent  test 
of  Doyle’s  expression  is  how  well  it  enables  the  coupled-dipole  method  to 
calculate  extinction  near  surface  mode  frequencies. 
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Table  A.l.  Comparison  of  Qext  for  nine  spheres  of  different  refractive 
indices  as  calculated  by  Mie  theory  and  by  the  coupled- 
dipole  method  with  either  the  Doyle  expression,  Draine’s 
radiation  reaction  term,  or  the  Clausius-Mosotti  relation. 
Values  shown  are  [Qext/Q„t(Mie)-l].  A  136-dipole  array  is 
used  for  the  coupled-dipole  calculations;  effective  radius 
of  the  sphere  remains  constant  at  3.19  du.  Rows  are  in 
order  of  increasing  size  parameter  of  the  dipolar  subunit. 


Refractive  Index 

Doyle 

Draine 

Clausius- 

Mosotti 

ka|m 

1.44  +  i0.26 

-1.0% 

-1.0% 

-1.0% 

0.29 

1.33  +  i0.05 

-1.7 

-3.3 

-5.0 

0.34 

1.55  +  10.005 

-6.8 

-11.1 

-12.8 

0.39 

1.39  +  i0.42 

0.0 

-1.1 

-0.6 

0.37 

1.7  +  iO.l 

-7.1 

-12.4 

-12.8 

0.42 

1.9  +  i0.0004 

-10.3 

-18.7 

-18.4 

0.48 

2.5  +  il.4 

14.1 

16.0 

16.7 

0.72 

3.5  +  i2.05 

17.3 

19.7 

15.3 

1.02 

3.0  +  i4.0 

30.0 

18.1 

36.9 

1.26 
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Surface  Modes 

Resonance  features  known  as  surface  modes  can  exist  in  small-particle 
absorption  spectra  even  when  none  occur  in  the  bulk  material.  Surface  modes 
are  caused  by  lattice  vibrations  in  insulators  and  most  commonly  occur  on  the 
high-frequency  side  of  bulk  absorption  bands.  The  absorption  characteristics 
associated  with  a  surface  mode  are  highly  shape  dependent  even  though  the 
particles  are  small  compared  with  the  wavelength. 

The  resonance  is  associated  with  the  vanishing  of  the  denominators  of 
the  Mie  scattering  coefficients.  In  the  limit  x  -*•  0  (finite  |m| ),  the 
denominator  of  an  vanishes  if  the  following  is  satisfied  (Bohren  and  Huffman, 
1983,  Sect  12.1): 


m 


2 


n+ 1 
n 


n  =  1,  2, .. 


(A.6) 


For  sufficiently  small  spheres  the  dominant  coefficient  is  at;  thus  for  n  =  1  the 
resonance  condition  is  m2  =  -2.  With  the  refractive  index  of  the  medium  close 
to  unity  the  complex  dielectric  function  of  the  particle  for  which  resonance 
occurs  is:  e  =  e'  +  ie"  =  -2  +  iO.  The  frequency  at  which  e'  =  -2  and 
e"  =  0  is  called  the  Frohlich  frequency;  the  corresponding  normal  mode  is 
known  as  a  Frohlich  mode.  Notice  also  from  equation  (1)  that  when  em  =  1, 
edip  =  -2  at  the  Frohlich  frequency.  Thus  according  to  the  CM  relation--either 
equation  A.l  or  A.3--the  polarizability  at  a  Frohlich  mode  is  infinite.  It  is 
under  this  condition  we  now  test  Doyle’s  expression  for  determining  the 
polarizability. 

An  example  of  a  Frohlich  mode  appears  in  the  infrared  for  quartz 
particles.  The  optical  constants  of  quartz  (Figure  A.4)  were  determined  by 
using  the  Lorentz  oscillator  model  and  published  dispersion  parameters 
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Figure  A.4  Complex  refractive  index  and  dielectric  function  as  a  function  of 
wavenumber  (cm'1)  for  quartz  as  obtained  from  the  Lorentz 
oscillator  model:  upper  curves  are  the  refractive  index, 

( - )  represents  the  real  part  and  ( - )  the  imaginary  part. 


Dielectric  Function 
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(Spitzer  and  Kleinman,  1961).  As  in  Huffman  and  Bohren  (1980)  we  used  the 
so-called  V3-%  rule  to  treat  the  anisotropy  of  quartz.  €'  attains  the  value  of  -2 
three  times--at  1153,  1165,  and  1169  cm^-with  attendant  values  of  e"  being 
0.39,  1.1,  and  0.49  cm'1,  respectively.  Although  e"  *  0  when  e'  =  -2,  e"  is 
sufficiently  small  to  produce  significant  maxima  in  the  extinction  cross  section, 
and  we  will  still  refer  to  these  as  Frohlich  modes.  As  seen  in  Figure  A.4,  the 
imaginary  part  of  the  refractive  index,  which  is  the  usual  first-order  indicator 
of  absorption,  undergoes  a  slight  inflection  near  1160  cm'1;  however,  the 
resonance  band  is  more  pronounced  than  this  indicates.  Having  all  the 
necessary  information  to  calculate  the  scattering  parameters  for  quartz 
particles,  our  next  step  was  to  compare  the  coupled-dipole  method  by  using  the 
Doyle  expression  and  the  CM  relation  one  last  time. 

Surface  Modes  for  Spherical  Particles 

Replacing  the  CM  relation  with  the  Doyle  expression  to  determine 
polarizability  in  the  coupled-dipole  method  leads  to  extinction  values  that 
compare  more  closely  with  Mie  theory.  This  was  shown  for  several  refractive 
indices  where  surface  modes  are  not  excited.  However,  when  the  calculations 
are  made  near  the  Frohlich  frequency,  the  CM  relation  does  better  than  the 
Doyle  expression.  The  following  model  comparison  shows  why. 

Extinction  efficiency  was  calculated  for  quartz  spheres  at  the  1153  cm'1 
Frohlich  frequency  by  using  Mie  theory  and  the  coupled-dipole  method.  The 
variable  of  interest  in  this  exercise  was  particle  radius,  which  ranged  from 
0.1  to  100  jum  when  modeled  by  Mie  theory.  The  coupled-dipole  method  was 
run  twice,  first  with  the  CM  relation  and  then  with  the  Doyle  expression,  for 
each  of  nine  spheres  with  radii  ranging  from  0.3  to  3.0  fim.  For  perspective, 
Mie  theory  was  also  used  to  calculate  Qext  for  two  other  refractive  indices. 

The  results  are  shown  in  Figure  A.5.  The  small-particle  resonance  associated 
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Radius  (pm) 


Figure  A3  Extinction  efficiency  as  a  function  of  particle  radius  (pim) 

obtained  from  Mie  theory  and  the  coupled-dipole  method.  Mie 
theory  is  represented  by  the  continuous  curves,  each  corresponds 

to  a  separate  complex  refractive  index:  - for  m  =  (0.1297, 

i 1.444)  (Frohlich  mode);  •  •  •  for  m  =  (133,  i0.05);  and - 

for  m  =  (1.7,  iO.l).  The  coupled-dipole  method  is  represented  by 
□  when  using  the  Clausius-Mosotti  relation  and  O  when  using 
Doyle’s  expression;  both  are  at  the  Frohlich  mode. 
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with  the  Frohlich  frequency  is  manifest  as  the  highest  Qext  in  the  radius  range 
from  3.0  /zm  and  smaller.  At  larger  radii  (>  3.0  /mi),  the  small-particle 
resonance  diminishes  and  becomes  obscured  by  the  bulk  properties  of  the 
sphere.  In  this  size  range,  values  of  Qext  for  all  three  refractive  indices 
approach  the  value  of  2.  It  is  in  the  size  range  between  0.3  and  3.0  /mi  for 
which  the  small-particle  resonance  is  pronounced,  and  Qext  obtained  by  using 
the  CM  relation  more  closely  agrees  with  Mie  theory  calculation  than  does  the 
Doyle  expression. 

It  is  observed  that  at  radii  of  0.3  and  3.0  /mi,  the  results  from  the 
coupled-dipole  method  with  the  CM  relation  or  the  Doyle  expression  are  in 
agreement.  It  has  been  stated  earlier  that  for  small  particles  the  Doyle 
expression  reduces  to  the  CM  relation,  which  is  why  the  two  coupled-dipole 
method  values  agree  at  0.3  /mi:  both  schemes  yield  identical  polarizabilities. 
For  particles  larger  than  3.0  /mi  the  bulk  extinction  features  become  more 
important  than  the  small-particle  resonance.  At  these  larger  particle  sizes  the 
CM  relation  and  the  Doyle  expression  give  different  values  for  polarizability, 
but  the  calculated  extinction  efficiencies  are  nearly  identical.  (Larger  spheres 
could  not  be  run  with  the  coupled-dipole  method  because  of  computer 
limitations.)  Thus  the  discrepancy  between  the  CM  relation  and  the  Doyle 
expression  lies  in  the  radius  range  from  0.3  to  3.0  /mi  and  appears  to  be  due  to 
the  small-particle  resonance.  This  is  not  unexpected  because  the  Doyle 
expression  contains  higher  order  terms  of  the  size  parameter  that  would  have 
their  greatest  effect  at  limiting  the  resonance  in  this  size  range  (according  to 
the  CM  relation,  polarizability  is  infinite  at  a  true  Frohlich  frequency,  but  it  is 
not  infinite  by  the  Doyle  expression).  Consequently  the  larger  polarizability 
from  the  CM  relation  leads  to  a  higher  Qext  as  seen  between  0.3  and  3.0  /mi. 

It  may  not  be  possible  for  the  coupled-dipole  method  to  duplicate  the 
resonance  peak  of  Mie  theory  at  the  Frohlich  frequency  because  the  dipolar 
array  cannot  exactly  represent  a  sphere.  (More  will  be  said  about  this  below.) 
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The  result  will  be  a  lower  resonance  peak.  Using  the  CM  relation  in  place  of 
the  Doyle  expression  yields  a  higher  polarizability  at  a  Frohlich  frequency,  so 
we  obtain  a  higher  resonance  peak.  Thus  conflicting  inadequacies  in  the 
model  permit  the  CM  relation  to  produce  extinction  values  at  the  Frohlich 
frequency  that  more  closely  agree  with  Mie  theory  than  does  the  Doyle 
expression. 

For  subsequent  modeling  the  Doyle  expression  will  be  used.  We  feel 
that  the  elimination  of  the  artificially  high  polarizability  obtained  with  the  CM 
relation  near  a  Frohlich  mode  should  be  avoided  despite  the  fact  that  the  CM 
relation  agrees  better  with  Mie  theory.  Replication  of  Mie  theory  is  a  goal 
when  modeling  spheres;  however,  in  reality  one  might  rather  model  a  so-called 
near-sphere  since  spheres  are  more  an  anomaly  in  nature  than  modeling 
practices  admit.  We  now  examine  extinction  by  quartz  particles. 

The  extinction  cross  section  per  unit  volume  Ccxt/v  for  a  spherical 
quartz  particle  of  radius  0.5  fxm  was  calculated  from  wavenumber  1000  to 
1300  cm'1  by  using  the  optical  constants  shown  in  Figure  A.4.  For  the  coupled- 
dipole  method,  an  array  of  136  dipolar  subunits  was  used  to  represent  the 
sphere.  The  results  appear  in  Figure  A.6.  Two  major  peaks  appear  in  the  Mie 
theory  results  that  correspond  to  the  Frohlich  frequencies--one  peak  is  at  1153 
cm’1  and  the  second  is  between  1165  and  1169  cm'1.  The  coupled-dipole 
method  also  predicts  extinction  peaks  near  these  frequencies;  however,  the 
peaks  are  not  as  high  as  those  computed  by  Mie  theory,  and  at  1153  cm1  the 
extinction  peak  is  shifted  toward  lower  frequency  and  has  a  broader  base. 

Both  issues  will  be  addressed;  first  is  the  extinction  magnitude. 

It  was  stated  earlier  that  the  resonance  at  a  Frohlich  frequency  is  highly 
shape  dependent.  How  well  can  a  collection  of  dipoles  on  a  cubic  lattice 
represent  a  sphere  at  this  frequency?  Draine  (1988)  presented  a  criterion  for 
estimating  the  sphericity  of  a  dipolar  array.  By  comparing  the  radius  of 
gyration  of  an  array  of  dipolar  subunits  with  a  sphere  of  equal  volume,  one 
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Wavenumber  (cm™1) 


Figure  A.6  Extinction  cross  section  per  unit  volume  (/mi'1)  as  a  function  of 
wavenumber  (cm'1)  for  spheres  as  obtained  from  Mie  theory 

( - )  and  from  the  coupled-dipole  method  using  the  Doyle 

expression  ( - ).  Peaks  correspond  to  the  Frohlich  frequency 

for  quartz. 


119 


can  pick  out  a  few  values  for  N  that  best  approximate  a  sphere.  We  adopted 
Draine’s  notation  of  [f(N)]w  as  the  ratio  of  the  radius  of  gyration  of  the  dipolar 
array  to  that  of  a  sphere.  Thus  as  a  measure  of  sphericity:  [f(N)]*  =  1  for  a 
perfect  sphere.  Draine  contended  that  surface  granularity,  which  is  inversely 
proportional  to  the  particle’s  sphericity,  leads  to  numerical  inaccuracies 
especially  when  the  refractive  index  is  large.  We  will  investigate  whether 
surface  irregularity  is  the  reason  the  coupled-dipole  method  cannot  duplicate 
the  peak  extinction  obtained  from  Mie  theory  at  the  Frohlich  frequency.  If 
this  were  the  case,  the  extinction  cross  section  could  be  affected  by  small 
changes  in  the  sphericity  of  a  dipolar  array. 

Using  the  radius  of  gyration  criterion,  Draine  considered  the  136-  and 
160-dipole  arrays  to  be  good  representations  of  spheres.  (He  also  considered 
several  larger  arrays  to  be  good  representations,  but  the  version  of  the 
coupled-dipole  used  here  is  limited  to  300  dipolar  subunits.)  These  dipole 
arrays  have  radii  of  3.0  du  and  3.5  du;  however,  the  central  dipoles  of  the 
arrays  are  displaced  from  the  origin  by  0.5  du.  With  the  central  dipole  located 
at  the  origin,  additional  spherical  approximations  are  obtained,  for  example: 
123-  and  179-dipole  arrays  (radii  are  3.0  and  3.5  du).  The  sphericity  of  seven 
arrays  was  computed  with  the  radius  of  gyration  criterion  and  is  listed  in 
Table  A.2.  Note  that  the  136-,  160-,  and  251-dipolar  arrays  have  radii  of 
gyration  closer  to  unity  than  the  remaining  four  arrays  which  implies  that  of 
the  seven  they  best  represent  a  spherical  particle. 

To  determine  if  small  changes  in  sphericity  affects  extinction  near  the 
1153  cm'1  Frohlich  frequency,  we  eliminated  the  size  difference  of  the  seven 
arrays.  The  radii  of  the  dipolar  subunits  of  the  seven  arrays  were  adjusted  so 
that  the  entire  array  represented  a  sphere  with  an  effective  radius  a,  of  0.5  pi m. 
Cext/v  was  then  calculated  for  the  seven  arrays;  the  values  are  reported  in 
Table  A.2.  At  non-surface  modes  Cext/v  for  spheres  is  unaffected  by  these 
variations  in  N  (when  ae  is  held  constant),  but  as  seen  in  Table  A.2  extinction 
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Table  A.2. 


Comparison  of  Ce„/v  among  seven  dipolar  arrays;  all 
represent  a  spherical  particle  of  radius  0.5  /*m.  Columns 
represent  number  of  dipoles,  radius  of  gyration,  and 
maximum  value  of  Cest/v  near  the  1153  cm'1  Frohlich 
frequency  as  calculated  by  the  coupled-dipole  method. 
Rows  are  in  order  of  decreasing  sphericity  as  determined 
by  Draine’s  criterion.  The  corresponding  CeJv  calculated 
by  Mie  theory  is  17.1  /mi'1. 


N 

[f(N)f 

Cew/V 

251 

1.0015 

14.40 

136 

1.0017 

13.61 

160 

1.0021 

14.94 

179 

1.0039 

10.93 

208 

1.0069 

11.58 

147 

1.0078 

9.44 

123 

1.0112 

10.12 
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values  differ  greatly  at  a  Frohlich  mode.  The  three  arrays  that  are  the  best 
representations  of  spheres  according  to  Draine’s  criterion  also  have  the  highest 
values  for  extinction.  This  is  not  meant  to  be  a  definitive  statement  but  it 
appears  that  extinction  is  dependent  on  how  well  the  dipolar  array  represents  a 
sphere.  Unfortunately,  larger  arrays  cannot  be  modeled  with  the  matrix 
inversion  method.  The  scattering  order  version  can  be  used  for  arrays  that 
contain  greater  than  1000  dipolar  subunits,  but  that  solution  method  diverges 
when  the  polarizability  is  too  high  (Singham  and  Bohren,  1988). 

A  corollary  of  this  is  that  bigger  arrays  may  not  always  be  better:  a 
sphere  should  not  necessarily  be  represented  by  the  largest  array  possible.  For 
example,  Table  A.2  implies  that  the  array  that  contains  136  dipolar  subunits 
represents  a  sphere  better  than  does  the  208-unit  array.  This  assertion  is  not 
conclusive,  and  needs  further  research. 

Although  sphericity  of  the  dipolar  array  affects  Cext/v,  it  does  not 
explain  why  values  are  lower  than  Mie  theory  by  as  much  as  20%.  This 
discrepancy  of  maximum  calculated  extinction  may  be  inherent  in  the 
formulation  of  the  coupled-dipole  method.  The  coupled-dipole  method  is  able 
to  match  the  extinction  computed  by  Mie  theory  at  the  Frohlich  frequency,  but 
only  when  the  dipolar  subunits  do  not  communicate  their  scattered  fields  to 
their  neighbors.  When  the  interactions  between  the  subunits  have  been  turned 
off,  they  scatter  as  individual  spheres  and  the  calculated  Cext/v  becomes 
identical  to  that  given  by  Mie  theory.  When  the  interactions  are  included,  the 
values  of  Cext/v  decrease  to  what  was  shown  in  Table  A.2.  Thus  the  coupling 
of  the  dipoles  appears  to  be  the  limiting  factor:  an  array  of  dipolar  subunits 
on  a  cubic  lattice  does  not  adequately  represent  a  sphere  when  the  shape 
dependence  is  as  crucial  as  at  a  Frohlich  mode.  We  now  will  investigate  the 
broadening  of  the  extinction  band  shown  in  Figure  A.6. 
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Shape  Effects  on  Surface  Modes 

Mie  theory  predicts  strong  resonance  peaks  at  Frohlich  frequencies  for 
quartz  crystals  in  the  Rayleigh  regime.  The  coupled-dipole  method  also  shows 
extinction  maxima  at  these  frequencies  as  noted  above;  however,  the  height 
and  width  of  the  coupled-dipole  extinction  curve  is  different  from  that 
obtained  with  Mie  theory.  The  reduction  in  height  was  explained  by  the 
inability  of  a  finite  array  of  dipolar  subunits  to  represent  a  sphere.  The 
widening  of  the  extinction  band  may  be  explained  in  a  similar  fashion. 

Resonance  features  are  not  limited  to  spheres.  Whereas  Mie  theory  is 
applicable  only  to  spheres,  the  coupled-dipole  method  may  be  used  to  calculate 
extinction  for  any  particle  shape.  To  investigate  nonspherical  resonances  near 
the  Frohlich  frequency,  we  modeled  the  following  shapes  by  using  the  coupled- 
dipole  method:  a  2x  1  and  a  4x  1  oblate  spheroid  and  a  2x  1  and  a  4x  1  prolate 
spheroid  comprised  of  152,  136,  138,  and  136  dipolar  subunits,  respectively. 

The  spheroids  were  modeled  with  an  equivalent  radius  of  0.5  fim  as  was  done 
when  modeling  the  136-dipolar  subunit  sphere  in  the  previous  section.  To 
simulate  random  orientation  of  the  particles  for  comparison  with  measured 
data  each  spheroid  was  rotated  in  over  100  alignments.  (The  spheres  in  the 
previous  section  were  also  rotated  but  little  effect  on  Qext  was  noted.) 

Values  of  Cext/v  for  the  four  spheroids  and  the  136-unit  sphere  are 
shown  in  Figure  A.7.  Near  the  Frohlich  frequencies,  the  spherical  array 
exhibits  the  highest  Cext/v.  Near  wavenumber  1125  cm'1,  a  larger  value  for 
extinction  is  calculated  for  both  oblate  spheroids  and  the  2  x  1  prolate  spheroid 
than  for  the  sphere  suggesting  the  presence  of  a  nonspherical  resonance. 
Recalling  Figure  A.6,  Cext/v  for  the  136-dipole  sphere  was  larger  at  1 125  cm'1 
than  that  computed  by  Mie  theory.  Since  the  136-dipole  sphere  was  shown 
earlier  to  be  slightly  nonspherical,  this  nonspherical  resonance  could  account 
for  this  larger  extinction  for  the  136-dipole  sphere.  The  proximity  of  the 
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Figure  A.7  Pictorial  representation  of  small-particle  extinction  for  quartz. 

Extinction  cross  section  per  unit  volume  calculated  by  the 

coupled-dipole  (Doyle  expression)  as  a  function  of  wavenumber 
(cm'1).  Five  shapes  are  represented:  O  symbolizes  a  sphere,  ||  a 
2x1  prolate  spheroid,  |  a  4x  1  prolate  spheroid,  =  a  2x  1 
oblate  spheroid,  and  —  a  4x  1  oblate  spheroid. 
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nonspherical  resonance  to  the  Frohlich  frequency  of  1153  cm'1  thus  widens  the 
extinction  band  of  the  136-dipole  sphere  towards  lower  wavenumbers  as  seen 
in  Figure  A.6.  A  similar  broadening  does  not  occur  towards  higher 
wavenumbers  because,  as  can  be  seen  from  Figure  A.7,  the  nonspherical 
resonances  at  1200  cm'1  are  weaker  than  those  at  1125  cm'1. 

The  results  from  modeling  the  sphere  and  spheroids  with  the  coupled- 
dipole  method  are  now  compared  with  another  theoretical  method--the  method 
of  continuous  distribution  of  ellipsoids  CDE  (Huffman  and  Bohren,  1980). 

This  method  is  used  to  compute  the  absorption  cross  section  per  unit  volume 
Cabs/v  by  integrating  over  a  distribution  of  shape  parameters  in  the  Rayleigh- 
ellipsoid  approximation.  The  CDE  expression  is  simplified  when  all  shapes 
within  the  ellipsoid  distribution  are  equally  probable: 


—  -  klmi-lS-Loge 
V  V  E  - 1  J 


(A.7) 


For  comparison  with  the  coupled-dipole  results,  extinction  is  taken  to  be  nearly 
equal  to  absorption.  This  is  not  a  bad  assumption  since  the  particles  are  small 
compared  with  the  wavelength,  and  in  the  vicinity  of  the  Frohlich  modes 
calculations  from  the  coupled-dipole  method  show  that  the  scattering  cross 
section  is  approximately  an  order  of  magnitude  smaller  than  the  absorption 
cross  section.  To  duplicate  the  equally-probable  spheroid  distribution  used  in 
the  CDE  method,  we  merely  averaged  the  Cext/v  values  from  the  coupled- 
dipole  method  for  the  sphere  and  spheroids.  The  modeling  results  are  shown 
in  Figure  A.8.  The  agreement  is  good  considering  the  CDE  method  represents 
a  continuous  distribution  of  spheroids  and  the  coupled-dipole  method  only  five 
arrays  (four  spheroids  and  a  sphere).  In  turn,  extinction  measurements  for 
submicron  quartz  crystals  have  been  shown  to  be  well  predicted  by  the  CDE 
method  (Huffman  and  Bohren,  1980).  Thus  by  the  transitive  property,  this 
exercise  has  given  credence  to  the  ability  of  the  coupled-dipole  method 
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Figure  A.8  Extinction  cross  section  per  unit  volume  (/mi’1)  as  a  function  of 
wavenumber  (cm'1)  obtained  from  the  Continuous  Distribution  of 

Ellipsoids  method  ( - )  and  the  coupled-dipole  method  (Doyle 

expression)  using  five  spheroids  (O):  sphere,  2x1  and  4x1  oblate 
and  prolate  spheroids. 
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with  the  Doyle  expression  to  treat  extinction  by  nonspherical  particles  near  a 
Frohlich  frequency. 

As  pointed  out  by  a  reviewer  the  coupled-dipole  method  may  be 
overkill  when  comparing  it  with  calculations  based  on  electrostatics.  However, 
because  retardation  is  accounted  for,  the  coupled-dipole  method  could 
(perhaps  should)  be  used  to  investigate  particles  with  larger  size  parameters 
near  a  Frohlich  mode  whereas  the  CDE  method  could  not.  The  purpose  of 
this  exercise  was  to  introduce  the  Doyle  expression.  As  for  comparing  the 
coupled-dipole  method  with  an  exact  method  such  as  the  T-matrix,  Goedecke 
and  O’Brien  (1988)  have  already  shown  agreement  between  the  two,  and, 
because  of  the  size  parameters  of  the  spheroid  as  particles  that  we  modeled, 
the  use  of  the  T-matrix  for  this  study  is  unwarranted. 

Conclusions 

The  coupled-dipole  method  relies  on  an  effective-medium  theory  to 
provide  the  polarizability  of  the  dipolar  subunits.  The  CM  relation  is  widely 
used  in  the  coupled-dipole  method,  but  because  of  the  size  of  the  dipolar 
subunits,  it  is  has  an  insufficient  number  of  terms  to  account  appropriately  for 
the  electric  dipole  polarizability.  Doyle  (1989)  used  the  electric  dipole  term 
from  Mie  theory  to  obtain  an  exact  expression  for  polarizability  that  afforded 
improved  calculations  of  optical  properties  of  a  suspension  of  metal  spheres. 
Incorporating  Doyle’s  expression  in  the  coupled-dipole  method  in  place  of  the 
CM  relation  yielded  improved  extinction  calculations  at  nonresonance 
frequencies. 

Doyle’s  expression  was  further  tested  by  calculating  extinction  by  quartz 
particles  near  a  Frohlich  mode.  The  small-particle  absorption  associated  with 
this  mode  is  highly  shape  dependent.  Mie  theory  and  the  coupled-dipole 
method  were  compared  by  calculating  the  extinction  by  spherical  particles; 
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differences  were  noted  and  examined.  First,  the  peak  extinction  calculated  by 
the  coupled-dipole  method  was  much  lower  than  that  obtained  from  Mie 
theory.  This  is  because  the  shape  dependency  at  the  Frohlich  frequency  is  so 
critical  that  a  spherical  array  of  dipoles  fails  to  represent  a  perfect  sphere. 
Second,  results  from  the  coupled-dipole  method  indicated  a  wider  absorption 
band  than  obtained  from  Mie  theory.  Further  modeling  with  the  coupled- 
dipole  method  indicated  the  broadening  arises  from  nonspherical  particle 
resonances.  The  calculated  extinction  by  these  nonspherical  particles  agrees 
with  the  results  of  another  theoretical  method,  the  CDE  method,  which  in  turn 
has  been  shown  to  agree  with  measured  data  (Hoffman  and  Bohren,  1980). 
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Appendix  B 


SUBROUTINES  FOR  BUILDING  DIPOLAR  ARRAYS 


This  appendix  contains  source  listings  of  the  subroutines  described  in 
Chapter  4.  The  subroutines  are  used  to  build  arrays  which  represent  particles  of 
several  different  shapes. 


non 
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SUBROUTINE  SPHERE  (RAD, REALRAD, OFFSET, X,Y,Z, ITOT, RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  A  SPHERE  OF  DIPOLAR  SUBUNITS. 

C 

C  INPUT:  RAD  IS  THE  NUMBER  OF  DIPOLAR  SUBUNITS  IN  THE  RADIUS  OF 
C  THE  SPHERE;  REALRAD  IS  THE  RADIUS  OF  THE- -IT  IS  USED  TO  ADJUST 

C  THE  SIZE  OF  THE  DIPOLAR  SUBUNITS;  OFFSET  -  0.0  IF  ORIGIN  IS 

C  LOCATED  WITH  THE  CENTRAL  DIPOLE,  OTHERWISE  -  0.5. 

C 

C  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATES  OF  THE  DIPOLAR 

C  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY;  RADEFF  IS 

C  THE  EFFECTIVE  RADIUS  OF  THE  SPHERE 
C 

DIMENSION  X(  *  ),  Y(  *  ),  Z(  *  ) 

OFFSET  =>0.0 

ITOT  =  0 

RAD2  -  RAD*RAD 

ISIZE  =  IFIX(RAD+1 . 0) 

C 

DO  30  JZ  -  -ISIZE, ISIZE, 1 
ZJZ  -  FLOAT (JZ)  +  OFFSET 
DO  20  JY  °  -ISIZE, ISIZE, 1 
YJY  -  FLOAT (JY)  +  OFFSET 
DO  10  JX  -  -ISIZE, ISIZE, 1 
XJX  =*  FLOAT  (JX)  +  OFFSET 

DETERMINE  IF  THIS  LATTICE  POSITION  FALLS  WITHIN  SPHERE'S  RADIUS 

R  -  1.0  *  (XJX**2  +  YJY** 2  +  ZJZ**2) 

IF  (R.GT.RAD2)  GO  TO  10 
C 

ITOT  -  ITOT  +  1 

X(ITOT)  -  1.0  *  XJX 

Y(ITOT)  -  1.0  *  YJY 

Z(ITOT)  =  1.0  *  ZJZ 

10  CONTINUE 

20  CONTINUE 
30  CONTINUE 
C 

C  RADEFF  =  ( (3*ITOT)/(4*PI ) )**(l/3)  IN  DIPOLE  UNITS 
RADEFF  =  (0.238732415  *  ITOT)**(l ./3 . ) 

C 

RAD  =  RADEFF 

IF  (REALRAD. NE. 0.0)  RAD  =  REALRAD 
C 

RETURN 

END 


o  n 
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SUBROUTINE  SPHROID  (AAA, BBB , REALRAD , X,Y, Z , ITOT .RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  AN  SPHEROID  OF  DIPOLAR  SUBUNITS. 

C 

C  INPUT:  AAA  AND  BBB  ARE  THE  DIMENSIONS  OF  THE  MAJOR  AND  MINOR  AXES; 

C  REALRAD  IS  THE  RADIUS  OF  THE  EQUIVALENT-VOLUME  SPHERE.  OFFSET  IS 

C  SET  AUTOMATICALLY.  CIRCULAR  CROSS  SECTION  IS  ORIENTED  IN  Y-Z 
C  PLANE.  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATES  OF  THE 

C  DIPOLAR  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY; 

C  RADEFF  IS  THE  EFFECTIVE  RADIUS  OF  THE  EQUIVALENT -VOLUME  SPHERE. 
C 

DIMENSION  X(  *  ),  Y(  *  ),  Z(  *  ) 

C 

IAA  -  IFIX(AAA  -  1.0) 

IBB  -  IFIX(BBB  -  1.0) 

DIAA  =  AAA/2. 

DIAC  -  BBB/2 . 

ITOT  -  0 
C 

DO  30  JX  -  -IBB, IBB, 2 
XJX  -  FLOAT (JX)/2 . 

DO  20  JY  -  -IAA, IAA, 2 
YJY  -  FLOAT( JY)/2 . 

DO  10  JZ  -  -IAA, IAA, 2 
ZJZ  -  FLOAT (JZ)/2. 

DETERMINE  IF  THIS  LATTICE  POSITION  FALLS  WITHIN  RADIUS 

RAD  -  ( (ZJZ/DIAA)**2  +  (YJY/DIAA)**2  +  (XJX/DIAC)**2) 

IF  (RAD. GT. 1.0)  GO  TO  10 
C 

ITOT  -  ITOT  +  1 
X(ITOT)  -  1.0  *  XJX 
Y(ITOT)  -  1.0  *  YJY 
Z(ITOT)  =  1.0  *  ZJZ 
10  CONTINUE 

20  CONTINUE 
30  CONTINUE 
C 

C  RADEFF  -  ((3*ITOT)/(4*PI))**(l/3) 

RADEFF  =  (0 . 238732415*ITOT)**(l ./3 . ) 

C 

RAD  -  REALRAD 

IF  (REALRAD. NE. 0.0)  RAD  =  REALRAD 
C 

RETURN 

END 


131 


SUBROUTINE  RECTSLD  (AAA, BBB , CCC .REALRAD ,X, Y, Z , ITOT .RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  A  RECTANGULAR  SOLID  OF  DIPOLAR  SUBUNITS. 

C 

C  INPUT:  AAA,  BBB,  CCC  ARE  THE  X,  Y,  Z,  DIMENSIONS;  REALRAD  IS  THE 
C  RADIUS  OF  THE  EQUIVALENT -VOLUME  SPHERE. 

C 

C  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATE  OF  THE  DIPOLAR 

C  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY;  RADEFF 

C  IS  THE  EFFECTIVE  RADIUS  OF  THE  EQUIVALENT -VOLUME  SPHERE. 

C 

DIMENSION  X(  *  ),  Y(  *  ),  Z(  *  ) 

C 

ITOT  -  0 

IAA  -  IFIX(AAA  -  1.0) 

IBB  -  IFIX(BBB  -  1.0) 

ICC  =  IFIX(CCC  -  1.0) 

C 

DO  30  JZ  =  -ICC, ICC, 2 
DO  20  JY  -  -IBB, IBB, 2 
DO  10  JX  -  -IAA, IAA, 2 
ITOT  -  ITOT  +  1 
Z(ITOT)  -  FLOAT (JZ)/2 . 

Y(ITOT)  -  FLOAT (JY)/2 . 

X(ITOT)  =  FLOAT (JX)/2. 

10  CONTINUE 

20  CONTINUE 
30  CONTINUE 
C 

C  RADEFF  -  ( (3*IT0T)/(4*PI) )**(l/3) 

RADEFF  -  (0.238732415  *  ITOT)**(l ./3 . ) 

C 

RAD  =  REALRAD 

IF  (REALRAD. NE. 0.0)  RAD  =  REALRAD 
C 

C  RETURN 

END 
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SUBROUTINE  CYLNDER  (RAD, HGHT .REALRAD , OFFSET, X, Y, Z , I TOT .RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  A  CYLINDER  OF  DIPOLAR  SUBUNITS. 

C 

C  INPUT:  RAD  IS  THE  NUMBER  OF  DIPOLAR  SUBUNITS  IN  THE  RADIUS  OF  THE 
C  CROSS  SECTION;  HGHT  IS  THE  LENGTH  OF  THE  CYLINDER;  REALRAD  IS  THE 

C  RADIUS  OF  THE  EQUIVALENT -AREA  CIRCLE  OF  THE  CROSS  SECTION; 

C  OFFSET  -  0.0  IF  THE  CENTRAL  DIPOLES  ARE  LOCATED  ON  THE  X  AXIS, 

C  OTHERWISE  =0.5.  CYLINDER  IS  ORIENTED  ALONG  X  AXIS. 

C 

C  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATES  OF  THE  DIPOLAR 

C  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY;  RADEFF  IS 

C  THE  EFFECTIVE  RADIUS  OF  THE  EQUIVALENT-AREA  CIRCLE  OF  CROSS 
C  SECTION. 

C 

DIMENSION  X(  *  ),  Y(  *  ),  Z(  *  ) 

ITOT  =  0 

RAD2  =  RAD  *  RAD 
ISIZE  =  IFIX(RAD  +  1.0) 

LSIZE  =  IFIX(HGHT  -  1.0) 

C 

DO  30  JY  =  -ISIZE, ISIZE, 1 
YJY  =  FLOAT (JY)  +  OFFSET 
DO  20  JZ  =  -ISIZE, ISIZE, 1 
ZJZ  =  FLOAT (JZ)  +  OFFSET 
R  =  (YJY**2  +  ZJZ**2) 

IF  (R.GT.RAD2)  GO  TO  20 
DO  10  JX  =  -LSIZE, LSIZE, 2 
ITOT  =  ITOT  +  1 
Z(ITOT)  =  ZJZ 
Y(ITOT)  =  YJY 
X(ITOT)  -  FLOAT( JX)/2 . 

10  CONTINUE 

20  CONTINUE 
30  CONTINUE 
C 

C  RADEFF=(ITOT/(HGHT*PI))**l/2  OF  THE  CYLINDER'S  RADIUS 

RADEFF  =  ( ITOT/(HGHT  *  3 . 141592654) )**. 5 
C 

C  RAD  =  RADEFF 

IF  (REALRAD. NE. 0.0)  RAD  =  REALRAD 
C 

RETURN 

END 
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SUBROUTINE  FCC  (RAD, REALRAD, OFFSET, X , Y, Z, ITOT, RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  A  SPHERE  OF  DIPOLAR  SUBUNITS  ON  A  FACE- 

C  CENTERED  CUBIC  LATTICE  RATHER  THAN  SIMPLE  CUBIC  LATTICE. 

C 

C  INPUT:  RAD  IS  THE  NUMBER  OF  DIPOLOAR  SUBUNITS  IN  THE  RADIUS  OF 

C  THE  SPHERE;  REALRAD  IS  THE  RADIUS  OF  THE  SPHERE;  OFFSET  =0.0 

C  IF  THE  ORIGIN  IS  LOCATED  WITH  THE  CENTRAL  DIPOLE,  OTHERWISE  =0.5. 
C 

C  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATES  OF  THE  DIPOLAR 

C  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY;  RADEFF  IS 

C  THE  EFFECTIVE  RADIUS  OF  THE  SPHERE. 

C 

DIMENSION  X(  *  ),  Y(  *  ),  Z(  *  ) 

ITOT  -  0 

ROOT2  =  SQRT ( 2 . ) 

FACE  -  ROOT 2  *  0.5 
RAD 2  -  RAD  *  RAD 
ISIZE  =  IFIX(RAD  +  2.0) 

C 

DO  50  JZ  =  -ISIZE, ISIZE, 1 
DO  40  JY  =  -ISIZE, ISIZE, 1 
DO  30  JX  -  -ISIZE, ISIZE, 1 

XJX  »  ROOT 2  *  FLOAT (JX)  +  OFFSET 
YJY  »  ROOT  2  *  FLOAT  (JY)  +  OFFSET 
ZJZ  =  ROOT2  *  FLOAT (JZ)  +  OFFSET 
DO  20  I  =  1,4 

R  -  1.0  *  (XJX**2  +  YJY**2  +  ZJZ**2) 

IF  (R.GT.RAD2)  GO  TO  10 
ITOT  =  ITOT  +  1 

X(ITOT)  =  1.0  *  XJX 

Y(ITOT)  =  1.0  *  YJY 

Z(ITOT)  -  1.0  *  ZJZ 

10  IF  (I.EQ.l)  THEN 

XJX  =  XJX  -  FACE 

YJY  -  YJY  -  FACE 

ELSE  IF  (I.EQ.2)  THEN 
XJX  =  XJX  +  FACE 

ZJZ  -  ZJZ  -  FACE 

ELSE 

XJX  =  XJX  -  FACE 

YJY  =  YJY  +  FACE 

END  IF 

20  CONTINUE 

30  CONTINUE 

40  CONTINUE 
50  CONTINUE 
C 

C  RADEFF  =  ((3*ITOT)/(4*PI))**(l/3) 

RADEFF  =  (0.238732415  *  ITOT/ROOT2)**(l ,/3 . ) 

C 

RAD  =  RADEFF 

IF  (REALRAD. NE. 0.0)  RAD  =  REALRAD 
C 

RETURN 

END 


non  no 
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SUBROUTINE  HEX  (AAA, BBB , REALRAD ,X, Y , Z , ITOT ,RADEFF) 

C 

C  THIS  SUBROUTINE  BUILDS  A  HEXAGON  SOLID  OF  DIPOLAR  SUBUNITS. 

C 

C  INPUT:  THE  PARTICLE  IS  ALIGNED  ALONG  THE  X  AXIS  WITH  TWO  OF  THE 
C  VERTICES  LYING  IN  THE  X-Z  PLANE.  AAA  IS  THE  NUMBER  OF  DIPOLES 
C  IN  THE  WIDTH  (Y  DIRECTION) ;  BBB  IS  THE  NUMBER  OF  DIPOLES  IN  THE 
C  LENGTH  (X  DIRECTION) ;  REALRAD  IS  THE  RADIUS  OF  THE  EQUIVALENT - 

C  AREA  CIRCLE  REPRESENTED  BY  THE  CROSS  SECTION  OF  THE  PARTICLE. 

C 

C  OUTPUT:  ARRAYS  X,  Y,  AND  Z  GIVE  THE  COORDINATES  OF  THE  DIPOLAR 

C  SUBUNITS;  ITOT  IS  THE  NUMBER  OF  DIPOLES  IN  THE  ARRAY;  RADEFF 

IS  THE  EFFECTIVE  RADIUS  OF  THE  CROSS  SECTION. 

DIMENSION  X(*),Y(*) ,Z(*) 

ITOT  -  0 

IAA  =  IFIX(AAA  -  1.0) 

IBB  -  IFIX(BBB  -  1.0) 

DETERMINE  Z  DISTANCE  WHICH  BEST  APPROXIMATES  EQUAL  LENGTH  FACES 

ICC  -  NINT(SQRT(2 . )  *  ( (AAA-1 . )/2  .) ) 

DO  30  JX=- IBB , IBB , 2 
C 

C  BUILD  THE  CENTER  SQUARE/RECTANGLE 
C 

IND  -  0 

DO  20  JY  =  -IAA, IAA, 2 
DO  10  JZ  -  -ICC, ICC, 2 
ITOT  -  ITOT  +  1 
Z(ITOT)  =  FLOAT (JZ)/2 . 

Y(ITOT)  -  FLOAT( JY)/2 . 

X(ITOT)  =  FLOAT (JX)/2 . 

10  CONTINUE 

20  CONTINUE 


BUILD  THE  UPPER  AND  LOWER  TRIANGLES 

DETERMINE  HOW  MANY  LINES  OF  DIPOLES  (NOL)  TO  ADD: 

NOL  -  IFIX(AAA/2 . ) 

DO  40  I  -  l.NOL 

START  BUILDING  NEW  LINES  REMEMBERING  TO  INDENT  ONE  EACH  TIME 

IND  -  IND  +  2 
JZ  -  ICC  +  IND 

DO  50  JY  *=  -  IAA+IND ,  IAA-  IND ,  2 
ITOT  -  ITOT  +  1 
Z(ITOT)  -  FLOAT (JZ)/2 . 

Y(ITOT)  -  FLOAT (JY)/2 . 

X(ITOT)  -  FLOAT (JX) /2 . 

ITOT  -  ITOT  +  1 
Z(ITOT)  -  FL0AT( JZ)/( -2 . ) 

Y(ITOT)  -  FLOAT (JY)/2 . 

X ( I TOT ) -FLOAT ( JX ) /2 . 

CONTINUE 

CONTINUE 

CONTINUE 

RADEFF=(ITOT/(BBB*PI))**l/2  OF  THE  CYLINDER'S  RADIUS 
RADEFF  =  (ITOT/(BBB*3. 141592654))**. 5 

RAD  =  RADEFF 

IF  (REALRAD .NE .0.0)  RAD  =  REALRAD 

RETURN 

END 
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